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I.  Summary 

Previous  reports  under  the  project  have  shown  that  in  order  to  make 
realistic  estimates  of  acoustic  gravity  wa"es  and  seismic  surface  waves 
generated  by  underground  events ,  such  as  explosions  and  earthquakes,  the  change 
in  surface  displacement  field  or  the  source  parameters  of  the  eve;^  must  be 
obtained.  Surface  displacement  data  for  actual  events  is  usually  limited 
and  of  dubious  quality. 

Section  II  deals  with  the  development  of  a  systematic  approach  of 
determining  the  dislocation  or  relative  displacement  on  the  fault  surface 
from  a  limited  set  of  observed  surface  displacement  data  with  their  esti¬ 
mated  errors.  The  technique  gives  not  only  a  best  fit  dislocation  model  but 
allows  one  to  determine  the  reliability  of  the  model  and  the  resolving  power 
of  the  data  set.  The  resulting  dislocation  model  can  then  be  used  to  extrap¬ 
olate  the  surface  displacements  outside  cf  the  data  set. 

-  I 

Section  III  deals  with  the  applications  of  the  inversion  technique 
described  in  the  previous  section  to  the  change  in  surface  displacement  data 
caused  by  deformation  resulting  from  the  1964  Alaska  earthquake.  For  this 
earthquake,  a  two-dimensional  finite  element  numerical  model  is  used  to 
calculate  surface  displacements  from  a  dislocation  imposed  on  a  fault  surface 
located  in  a  heterogeneous  medium.  The  inversion  technique  is  used  to 
calculate  a  dislocation  model  which  fits  the  observed  data  to  a  high  degree 
of  accuracy.  An  error  analysis  is  carried  out  for  the  plain-strain  approxi¬ 
mation,  and  the  resolvability  of  the  features  of  the  calculated  dislocation 
is  examined.  The  results  indicate  that  the  observed  deformation  occurred 


as  the  result  of  massive  underthrusting  of  the  Alaska  continental  block  by 
the  dovmgoing  Pacific  plate. 

Sections  II  and  III  are  chapters  II  and  III  of  Ralph  Wilson  Alewine,  Ill's 
Thesis  (1974)  titled  "Application  of  Linear  Inversion  Theory  Toward  the 
Estimation  of  Seismic  Source  Parameters".  References  and  details  cited  in 
these  sections  can  be  found  in  the  Thesis.  This  thesis  has  been  submitted 
and  successfully  defenued  as  partial  fulfillment  of  the  requirements  for  the 
Degree  of  Doctor  of  Philosophy  at  the  California  Institute  of  Technology, 


Pasadena,  California. 
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Chapter  2 

Development  of  the  Stochastic  Inverse 
as  Applied  to  Static  Dislocation  Problems 

2.1  Introd uctlon. 

The  first  requisite  in  the  application  of  any  theory 
toward  the  estimation  of  seismic  source  parameters  is  the 
ability  to  solve  the  forward  problem  for  the  observed 
data  type  for  an  appropriate  seismic  sourct.  This  means 
simply  that  given  a  certain  method  of  physically  describ¬ 
ing  a  source  (analytically,  numerically,  or  by  analogue) 
we  are  able  to  estimate  changes  in  data  for  a  given  value, 
or  change  in  value,  of  particular  parameters  which 
describe  the  source.  Mathematically  this  is  mapping 
changes  in  the  source  model  space  into  changes  in  the 
data  space.  What  will  be  discussed  first  in  this  chapter 
is  Just  this  process,  and  later  we  will  look  at  the 
inverse  of  this  process.  By  the  inverse  of  this  process, 
we  mean  that  given  some  observations,  what  estimates  can 
be  made  about  the  different  source  parameters  which 
describe  our  source? 

The  Inversion  scheme  for  static  data  that  we  propose 
in  this  chapter  has  the  provision  for  the  inclusion  of 
the  estimated  variance  of  the  data  that  is  to  be  inverted. 
The  inclusion  of  this  data  variance  gives  rise  to  the 
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fact  that  we  cannot  now  ■‘Stimate  the  seismic  source  para¬ 
meters  exactly,  since  small  perturbations  in  the  model 
source  parameters  might  cause  changes  in  the  calculated 
data  values  which  11'..  inside  the  estimated  data  error 
limits.  This  concept  gives  rise  to  our  wondering  what 
ability  that  we  have  to  actually  resolve  any  detail  of  the 
various  parameters  of  our  fault  model.  This  resolution 
question  will  be  examined  in  some  detail  in  this  chapter. 

Vve  will  first  consider  the  problem  of  estimating 
source  parameters  for  static  data.  The  procedure  devel¬ 
oped  for  this  case  can  then  be  extended  to  t.iat  of  esti¬ 
mating  dynamical  source  parameters.  This  extension  is 
done  in  a  later  chapter.  A  brief  review  of  the  develop¬ 
ment  of  static  field  solutions  due  to  various  earthquake 
sources  is  in  order. 

2 .2  Development  of  the  Forward  Static  Problem . 

Numerous  attempts  have  been  made  in  the  past  several 
years  to  interpret  the  observed  permanent  changes  in  the 
displacement  and  suraln  fields  due  to  the  occurrence  of  an 
earthquake.  Various  approaches  to  the  solution  of  this 
problem  have  been  proposed,  each  based  on  a  slightly 
different  interpretation  of  the  earthquake  source  process 
as  a  whole.  In  each  of  these  approaches,  there  exists  in 
the  interior  of  the  elastic  medium  some  discontinuity 
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surface  which  we  can  associate  with  a  "fault".  The  diff¬ 
erent  initial  or  boundary  conditions  that  can  be  applied 
to  this  discontinuity  surface  give  rise  to  the  various 
approaches.  These  approaches  can  be  broken  into  four  main 
groups:  stress  pulse  theory,  stress  relaxation  theory, 
dislocation  theory,  and  numerical  analogues.  With  the 
exception  of  the  dislocation  theory,  which  wil  be  treated 
in  more  detail,  a  short  description  of  the  approach  of 
each  of  these  theories  will  be  presented.  The  dislocation 
theory  is  reviewed  in  more  detail  because  it  involves  a 
parameter  that  is  readily  observable  when  the  discontin¬ 
uity  surface  breaks  the  free  surface  —  a  physical  offset. 
In  addition,  it  is  somewhat  more  straightforwardly  pleas¬ 
ing  to  model  static  dislocations  on  the  surface  caused  by 
static  dislocations  imposed  within  the  medium  rather  than 
the  more  obscure  parameter j  —  stress  and  strain. 

However,  it  will  be  seen  that  the  other  theoretical 
approaches  can  be  equivalenced  to  so™°  dislocation  repre¬ 
sentation  in  the  static  limit. 

Dislocation  Theory.  As  a  mathematical  mod.el  of  a  "fault", 
the  concept  and  formulation  of  a  physical  dislocation  has 
been  extensively  used.  The  dislocation  surface  in  an 
elastic  medium  is  viewed  as  a  surface  over  which  there  is 
a  discontinuity  in  displacement.  One  of  the  first  efforts 
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to  explain  the  elastic  displacements  resulting  from  a 
dislocation  was  that  done  by  Vvendenskaya  (1956,  1958). 
Probably  the  most  lucid  explanation  of  the  dislocation 
theory  for  calculating  static  changes  that  accompany 
faulting  was  given  in  a  set  of  papers  by  Steketee  (1958a, 
1958b).  In  these  papers,  Steketee  recognized  that  the 
relations  for  the  displacement  field  in  an  infinite  elas¬ 
tic  medium  strained  by  a  dislocation  over  some  surface  as 
given  earlier  Dy  Volterra  (1907)  would  be  appropriate  in 
describing  the  deformation  that  accompanies  faulting. 
Steketee  derived.,  through  the  use  of  Galerkin  vectors,  the 
expressions  for  static  displacements  in  an  infinite 
elastic  solid.  These  relations  were  given  in  compact  form 
as  integrals  over  the  dislocation  surface. 

In  his  papers,  Steketee  poses  the  following  problem. 

A  dislocation  surface,  X  ,  is  created  with  an  elastic 
solid  which  is  bounded  by  some  surface  S  .  The  media  is 
then  strained  by  the  introduction  of  a  certain  distribu¬ 
tion  of  "nuclei  of  strain"  (Love,  19M)  along  the  dislo¬ 
cation  surface.  The  nuclei  were  shown  to  exist  in  si:: 
basic  forms  corresponding  broadly  to  a  combination  of  a 
center  of  dilatation  and  a  double  force  without  moment, 
and  secondly,  two  coplanar,  mutually  perpendicular  double 
forces  with  moments.  For  pure  shear  dislocations,  only 


the  latter  type  nuclei  are  applicable,  whereas,  for  pure 
tensile  dislocations,  the  former  are  applicable.  The 
displacements  at  a  point  Q,  u^(Q) ,  within  the  elastic 
solid  can  be  written  as 


uk«S) 


sW/v  iui(p)“ij(p-Q,vjd2 


*  mffs  ui“ij(p’9)vjds 


(2.1) 


In  this  equation  v4  are  the  direction  cosines  of  the 

o 

normal  to  the  dislocation  surface  elements,  y  is  the 

rigidity,  and  Au^(P)  is  the  dislocation  function  on  the 

surface  2.  It  is  seen  that  for  an  arbitrary  dislocation, 

a  set  of  six  of  these  functions  is  necessary.  (i=l,2,3 

j=l,2,3  with  ij=ji.)  The  kernels  of  the  integrals, 
k 

w^j(P,Q),  are  the  displacements  at  the  observation  point 
due  to  a  single  nucleus  of  strain  A  summation  over  all 
nuclei  is  implied.  As  the  surface  S  is  enlarged  to  in¬ 
finity,  the  displacements,  u^ ,  on  S  approach  zero  and  the 
second  integral  vanishes. 

The  formalism  for  this  problem  was  extended  to  include 


a  dislocation  in  a  semi-infinite  elastic  medium  by  a 
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superposition  of  solutions,  which  together  satisfy  the 
required  boundary  conditions  at  the  free  surface.  These 
boundary  conditions  require  the  solution  to  be  a  fairly 
complex  boundary  value  problem,  however,  it  is  cleverly 
solved  by  a  superposition  of  solutions  in  the  following 
manner.  The  tangential  stress  at  the  free  surface  is 
made  to  vanish  by  the  addition  of  an  image  dislocation 
"above”  the  free  surface.  This  last  superposition  is 
commonly  referred  to  as  the  Boussinesq  problem.  The 
strength  of  the  Boussinesq  load  is  such  to  cancel  the 
normal  stress  on  the  free  surface  which  is  doubled  by  the 
addition  of  the  image  source.  Using  the  Volterra  rela¬ 
tions,  the  displacement  field  at  a  point  Q  in  a  semi- 
infinite  medium  is  then  given  by 

VQ>  *  *u1(P)^J(P.«l)»J«  •  (2-2) 


Comparison  of  (2.1)  and  (2.2)  shows  that  only  the  values 
of  the  kernels  are  changed  by  the  imposed  boundary 
conditions.  The  kernels  of  (2.2),  ,  are  the  set  of 

Green's  functions  found  from  the  superposition  of  solu¬ 
tions  which  satisfy  these  half-space  boundary  conditions. 
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Steketee  (1958a)  gives  the  exact  form  of  one  of  these 
functions,  W*2,  which  is  all  that  is  necessary  to  approxi¬ 
mate  a  vertical  strike  slip  fault.  Chinnery  (1961,  1963, 
1965)  took  the  general  expression  (2.2)  and  derived  an 
exact  analytical  form  of  the  displacement  and  stress  fields 
on  the  surface  of  a  semi— infinite  medium  for  an  internal 
rectangularly-shaped  dislocation  surface  modeling  a  verti¬ 
cal  strike  slip  (transcurrent)  fault.  In  performing  these 
calculations,  Chinnery  assumed  that  the  dislocation  dis¬ 
continuity  was  constant  over  the  entire  fault,  and  he  also 
assumed  that  the  Lame  parameters  for  the  solid  were  equal 
so  that  the  integration  could  be  carried  out  exactly. 

Thus,  the  elastic  medium  for  which  this  theory  is  applica¬ 
ble  is  one  in  which  the  Poisson  ratio  is  constant  at  0.25. 
Steketee  (1958a)  showed,  however,  that  (2.2)  is  valid  where 
Au^(P)  takes  any  form  (Somigliana  dislocation)  provided 
that  the  tensile  forces  across  the  dislocation  surface  sum 
to  zero. 

Maruyama  (1964)  has  derived  the  remaining  five  sets 
of  Green’s  functions  needed  in  the  solution  of  an  arbi¬ 
trary  dislocation  problem.  He  further  gives  explicit , 
analytic  solutions  for  the  displacement  field  at  the  free 
surface  due  to  constant  finite  dislocations  on  rectangular 
surfaces.  The  dislocations  considered  are  those  only  along 
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single  primary  axes.  The  examples  th\t  he  presents  in¬ 
clude  cases  for  which  the  dislocation  surface,  I,  is  both 
perpendicular  and  parallel  to  the  free  surface. 

Maruyama  (1963)  and  Burridge  and  Knopoff  (1964) 
showed  that  the  displacement  fields  produced  by  a  disloca¬ 
tion  on  a  mathematical  description  of  a  dislocation  fault 
surface  is  equivalent  to  that  produced  by  a  suitable  dis¬ 
tribution  of  forces  on  the  fault  surface  acting  as  if  there 
was  no  fault  present.  Utilization  of  this  fact  makes 
possible  the  use  of  work  in  mathematical  elasticity  theory 
done  much  earlier  than  Steketee's  (1958a)  work.  Notable 
in  this  early  literature  is  that  by  Mindlin  (1936)  who 
treated  the  static  problem  of  a  single  force  acting  in  a 
homogenous  half-space.  Mindlin  and  Cheng  (1950)  give 
explicit  expressions  for  the  displacement  and  stress  fields 
due  to  point  forces  and  double  forces  acting  in  an  elastic 
half-space.  Maruyama  (1964)  gives  a  short  summary  of  the 
early  literature  in  Japan  and  elsewhere  on  this  subject. 
This  include"!  work  done  by  Sezawa  (1929)*  Honda  and  Miura 
(1935),  Whipple  (1936),  SoeJa  (1944)  and  Yamakawa  (1955)- 
Press  (1965)  showed  that  the  kernels  of  (2.2)  could  be 
derived  in  a  straightforward  manner  from  the  results  of 
Mindlin  and  Cheng  (1950).  Press  obtained  the  same  results 
for  a  vertical  strike  slip  fault  as  Chinnery  had  lone 
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prevlously,  and  obtained  the  same  results  for  a  vertical 
dip  slip  fault  as  Maruyama  (1964)  derived.  In  this  paper. 
Press  (1965)  addeu  the  analytic  expressions  for  tilts  and 
strains  for  these  particular  fault  orientations.  Savage 
and  Has  tie  (1966)  used  the  theory  given  by  Maruyama  (1964) 
to  calculate  the  vertical  displacements  induced  by  dis¬ 
locations  on  fault  surfaces  that  could  have  components  of 
dip  other  than  in  a  direction  perpendicular  to  the  free 
surface.  This  led  to  the  ability  to  model  more  geologi¬ 
cally  realistic  faults. 

Mansinha  and  Smylie  (1971)  completed  the  derivation 
of  the  displacement  fields  due  to  buried  dislocations  on 
finite  rectangular  surfaces.  These  authors  give  the  com¬ 
plete  closed  form,  indefinite  integral  expressions  for 
the  entire  displacement  fields,  both  at  the  free  surface 
and  at  any  depth  in  the  elastic  half-space,  due  to  a 
rectangular  dislocation  surface  that  can  be  arbitrarily 
inclined.  The  fields  are  presented  in  such  a  form  that 
they  are  readily  evaluated  numerically  on  the  computer  and 
Involve  only  simple  algebraic  and  trigonometric  functions. 
However,  these  authors  do  not  give  the  formulas  for  the 
strain  and  tilt  fields  arising  from  a  dislocation  across 
an  arbitrarily  inclined  surface.  These  strain  and  tilt 
fields  can  be  easily  obtained  from  differentiation  of  the 
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displacement  fields.  Appendix  1  of  this  thesis  gives  the 
results  of  this  differentiation. 

Chinnery  and  Petrak  (1968)  extended  the  work  of 
Chinnery  (1961)  by  considering  a  model  of  a  vertical 
strike  slip  fault  on  which  the  dislocation  uniformly  and 
smoothly  goes  to  zero  near  the  edges  of  the  dislocation 
surface.  This  variation  was  chosen  so  as  to  remove  the 
stress  singularity  that  was  occurring  at  the  edge  of  the 
fault  surface  in  the  original  work.  Except  in  extreme 
cases,  the  tapering  of  the  dislocation  near  the  edges  of 
the  surface  had  little  effect  on  the  overall  displacement 
fields  calculated  on  the  surface. 

Ben-Menahem  and  Singh  (1963a)  and  Ben-Menahem  and 
Gillon  (1970)  computed  the  integral  expressions  for  the 
displacement  field,  both  dynamic  and  static,  at  the  free 
surface  for  a  model  of  a  vertical  strike  slip  fault  and  a 
vertical  dip  slip  fault  for  a  medium  which  contains  a 
layer  of  arbitrary  thickness  over  a  uniform  half-space. 
These  authors  point  out  that  due  to  the  complexity  of  the 
problem,  the  use  of  the  Galerkin  vectors  for  elastic 
problems  involving  more  than  one  layer  over  a  half-space 
would  be  extremely  difficult.  These  authors  suggested  the 
use  of  a  method  employing  Hansen’s  eigenvectors  in  obtain¬ 
ing  the  static  response  of  a  multilayered  homogenous 
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half-space.  McGinley  (1969)  and  Sato  (1971)  achieved 
much  the  same  results  as  these  authors  by  the  superposi¬ 
tion  of  several  half-space  Green’s  functions  solutions 
off-set  in  such  a  way  as  to  represent  a  layered  half¬ 
space.  Braslau  and  Lieber  (1968)  solved  the  static 
linearly  elastic  problem  of  a  concentrated  vertical 
Volttrra  dislocation  in  a  layer  over  a  half-space  .  They 
made  use  of  a  special  displacement  function  which  they 
called  a  modified  Galerkin  vector  to  give  the  solution  in 
a  form  which  must  be  evaluated  numerically.  Singh  (1970, 
1971)  has  applied  the  Thomson-Haskell  matrix  propagation 
method  (Thomson,  1950;  Haskell,  1953)  to  solve  the  problem 
of  static  deformation  in  a  multilayered  elastic  half-space. 
He  obtains  source  functions  for  the  six  elementary  dis¬ 
locations  that  were  given  by  Steketee  (1958a) .  Explicit 
integral  expressions  are  given  for  the  surface  displace¬ 
ments  for  a  vertical  strike  slip  and  vertical  dip  slip 
fault  when  these  faults  can  be  represented  by  concentrated 
or  point  sources.  Extension  to  finite  size  sources  is 
given  as  another  Integration  involving  the  dislocation 
surface.  Recently  Chinnery  and  Jovanovich  (1972)  have 
calculated  the  displacement  field  due  to  a  vertical 
strike-slip  fault  of  infinite  length  for  an  earth  model 
consisting  of  two  layers  of  arbitrary  thickness  and 
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rigidity  over  a  half-space.  Their  expressions  are  given 
in  series  form  so  that  no  further  integration  is  necessary. 
On  the  basis  of  this  model,  they  conclude  (and  thus  agree 
with  McGinley  (1969))  that  the  presence  of  a  low  rigidity 
layer  would  have  a  very  strong  (amplifying)  effect  on  thj 
observed  displacements  and  strains  in  the  far  field. 

Ben-Menahem  and  Singh  (1968b)  treated  in  detail  the 
problem  of  deformation  of  a  uniform  non-gravitating  sphere 
due  to  internal  Volterra  type  dislocations  of  arbitrary 
orientation  and  depth.  This  work  was  subsequently  expanded 
(Ben-Menahem  et_  al. ,  1969,  1970;  Singh  and  Ben-Menahem, 
1969;  Ben-Menahem  and  Singh,  1970;  Wason  and  Singh,  1972) 
to  include  the  computations  for  the  displacement  and 
strain  fields  everywhere  on  the  surface  of  a  homogenous 
sphere  induced  by  an  internal  dipolar  source  of  finite 
size.  The  results  for  a  sphere  were  shown  to  be  quite 
different  than  that  expected  in  the  far-field  half-space 
problem. 

Stress  Pulse  Theory.  This  approach  has  seen  limited  use 
in  explaining  elastostatic  phenomenon.  Kasahara  (1957) 
devised  this  method  to  model  the  mechanism  of  an  earth¬ 
quake  as  a  distribution  of  stresses  or  strains  imposed  on 
an  underground  plane.  When  the  conditions  of  elastic 
equilibrium  are  satisfied,  the  deformations  at  the  surface 
can  be  calculated.  He  models  an  infinite  strike-slip 


fault  with  a  zone  of  constant  stress  extending  to  a  given 
depth.  The  faulting  occurs  by  the  liberation  of  this 
initially  applied  shear  stress.  Horizontal  displacements 
were  calculated  for  various  depth  extensions  and  compari¬ 
sons  were  made  to  actual  faults  .  By  examining  the  diminu¬ 
tion  of  horizontal  displacement  with  distance,  the  depth 
of  extension  of  this  constant  stress  zone  is  determined. 
This  mechanism  is  extended  in  a  second  paper  (Kasahara, 
1959)  to  include  non-vertical  strike-slip  faults.  The 
static  mechanism  presented  by  Kasahara  is  analogous  to 
the  stress  pulse  problems  encountered  in  dynamical  formu¬ 
lations  of  seismic  sources.  Minster  (197*0  describes  the 
mathematical  nuances  of  this  approach. 

Stress  Relaxation  Theory.  A  third  method  of  determining 
the  static  deformation  from  a  model  of  an  earthquake  is 
obtained  through  an  entirely  different  approach  to  the 
theoretical  problem.  The  methods  considered  thus  far  are 
all  based  on  relations  in  which  conditions  on  various 
boundaries  are  imposed  (boundary -value  problems). 
Archambeau  (1964,  1968)  has  proposed  sin  alternative  mech¬ 
anism  of  describing  the  processes  which  accompany  the 
occurrence  of  earthquakes  —  that  of  material  failure. 
This  theory  is  devised  in  the  context  of  an  initial-value 
problem  in  that  a  medium  is  assumed  to  be  initially  in 
some  prestressed  state.  Deformation  in  the  medium  is 
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caused  by  introducing  some  surface,  or  volume,  within  the 
medium  where  the  material  falls.  This  failure  is  accomp¬ 
lished  by  making  a  significant  reduction  in  the  shear 
tractions  across  the  failure  surface.  The  media  then 
responds  by  "relaxing”  to  a  new  equilibrium  state  by 
radiating  the  energy  released  from  the  local  reduction  in 
strain  energy  in  the  source  region.  This  theory  has  been 
very  successful  in  the  dynamical  regime,  most  notably  in 
the  prediction  of  far-field  radiation  patterns  from  earth¬ 
quakes  and  explosions  accompanied  by  tectonic  release 
(Archambeau  and  Sammls ,  1970;  Lambert  et  al . ,  1972; 
Archambeau,  1972).  Because  of  the  theoretical  complexi¬ 
ties,  this  source  formulation  has  not  yet  been  directly 
applied  to  near-field  static  deformation  problems. 

Minster  (197*0  has  discussed  from  a  mathematical 
point  of  view  in  some  detail  the  similarities  and  differ¬ 
ences  between  the  various  formulations  of  the  earthquake 
processes.  Although  his  approach  is  mainly  based  on 
dynamical  considerations,  he  shows  that  in  the  static 
limit  the  general  representation  of  the  stress  relaxation 
and  stress  pulse  problems  reduce  to  the  displacement  field 
as  given  by  a  generalized  Somigliana  dislocation  along  a 
surface  of  shear  displacement  discontinuity.  This  same 
proof  was  attempted  by  McGinley  (1969),  but  the  arguments 
presented  by  Minster  (197*0  are  much  more  compxete. 
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Therefore,  we  may  express  the  source  in  terms  of  a 
Somigliana  dislocation  without  loss  of  equivalency  from 
the  other  source  descriptions.  An  approximation  to  this 
Somigliana  source  will  be  adopted  throughout  this  thesis. 
Numerical  Analogue.  An  altogether  different  approach  to 
solving  the  forward  problem  for  dislocations  in  an  elastic 
half-space  is  afforded  through  the  use  of  the  finite 
element  numerical  technique.  Use  of  this  technique,  which 
usually  requires  a  large  computing  capability,  enables 
solutions  to  be  found  to  problems  involving  heterogenities , 
both  later  and  vertical,  and  anisotropy  just  as  easily  as 
those  involving  a  uniform  homogenous,  isotropic  half-space. 
The  mechanics  of  this  method  have  been  described  exten¬ 
sively  in  the  engineering  literature  (Martin,  1966; 
Przemienlckl ,  1968;  ^.enkins,  1969;  Zienkiewicz,  1971). 

In  this  technique,  the  elastic  half-space  continuum  is 
divided  into  geometric  elements  which  are  inter-connected 
only  at  a  finite  number  of  nodal  points.  It  is  at  these 
nodal  points  that  displacements,  stresses,  or  forces  can 
be  imposed  on  the  system.  Concurrently,  stresses  and 
displacements  at  a  distance  removed  from  these  disturban¬ 
ces  can  only  be  measured  at  these  nodal  points.  The 
solution  to  the  system  of  simultaneous  equations  generated 
by  a  disturbance  imposed  on  a  given  node  is  constrained 
by  the  boundary  conditions  relevant  to  the  problem  and  is 
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solved  numerically.  Jungels  (1973)  gives  a  description  of 
the  adaptation  of  this  method  to  the  modeling  of  disloca¬ 
tion  fault  surfaces.  The  reader  is  referred  to  this  work 
for  a  summary  of  the  intricacies  of  this  numerical  method. 

Jungels  (1973)  and  Jungels  and  Frazier  (1973)  make  a 
positive  comparison  between  the  calculated  static  dis¬ 
placement  field  due  to  a  dislocation  in  a  uniform  homoge¬ 
nous  elastic  half-space  calculated  by  the  finite  element 
method  and  by  the  conventional  exact  Green's  functions 
techniques.  Although  this  author  had  at  his  disposal  a 
numerical  code  which  would  allow  only  the  modeling  of 
plane  strain  problems,  i.e.,  faults  of  infinite  length, 
more  recent  finite  element  numerical  codes  can  accommodate 
problems  involving  finite  dimensions  in  all  directions. 

The  great  advantage  of  this  method  in  calculating  dis¬ 
placement  and  strain  fields  from  models  of  earthquakes  is 
the  ability  to  vary  the  elastic  properties  of  the  medium 
both  over  the  fault  surface  and  the  source  to  observer 
path.  This  technique  can  be  limited,  however,  by  the 
shear  size  of  computer  storage  necessary  to  solve  a  prob¬ 
lem  in  which  the  continuum  must  be  very  finely  sampled 
in  order  to  accurately  approximate  the  continuum  for  the 
order  of  the  disturbance  being  modeled. 
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2.3  Application  of  the  Forward  Problem  coward  the 

Explanation  of  Observed  Static  Data. 

As  Is  obvious  from  the  preceding  discussion,  much 
progress  has  been  made  toward  the  static  modeling  of  the 
earthquake  source.  The  state-of-the-art  is  such  that  now 
an  accurate  description  of  the  static  processes  accompany¬ 
ing  faulting  can  be  investigated.  However,  tne  inverse 
problem  now  remains.  As  the  facility  for  calculating  the 
displacement  and  strain  fields  from  fault  models  became 
more  sophisticated,  a  wider  range  of  cat  a  came  under 
scrutiny  in  trying  to  infer  some  information  about  the 
various  parameters  which  affect  the  faulting  process.  The 
earliest  attempt  to  extract  source  information  from  static 
data  was  applied  to  differential  horizontal  displacements 
measured  near  long  vertical  strike-slip  faults.  Kasahara 
(1957,  1959),  Chinnery  (1961),  and  Chinnery  and  Petrak 
(1968)  tried  to  infer  the  depth  and  distribution  with 
depth  of  dislocation  faulting  by  fitting  the  rate  of  fall- 
off  of  horizontal  displacements  measured  parallel  to  the 
fault  strike  as  a  function  of  distance  away  from  the  sur¬ 
face  expression  of  the  fault.  A  trial  and  error  method 
was  used  to  fit  the  data  and  to  try  to  exclude  possible 
faulting  models.  Press  (1965)  and  Press  and  Jackson  (1965) 
used  Press'  calculations  to  model  the  close-in  vertical 
movements  associated  with  the  196*4  Alaskan  earthquake. 
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These  movements  were  modeled  w^th  a  vertical  dip-slip 
fault  and  an  indication  was  sought  as  to  the  depth  of 
faulting.  A  comparison  of  the  calculated  and  observed  far- 
field  residual  strain  steps  was  also  undertaken.  Singh  and 
Ben-Menahem  (1969)  attempted  to  fit  the  same  strain  obser¬ 
vations  using  their  meth  for  taking  into  account  the 
earth's  curvature.  In  both  these  studies,  no  attempt  was 
made  to  systematically  vary  the  source  parameters  to 
achieve  the  best  fit  to  the  data. 

As  displacement  data  foi  large  earthquakes  became  more 
abundant  and  reliable,  it  became  apparent  that  the  simple 
fault  models  having  a  constant  dislocation  over  the  entire 
fault  surface  could  not  adequately  represent  the  observa¬ 
tions.  Stauder  and  Bollinger  (1966)  first  proposed  that 
differential  slip  on  the  fault  surface  might  provide  a 
more  realistic  model  to  better  fit  the  data  from  the  1964 
Alaskan  earthquake.  They  approximated  the  differential 
movement  by  allowing  the  displacement  on  the  fault,  Au,  to 
vary  piecewise  along  the  direction  of  the  slip.  To  do 
this,  the  total  fault  plane  was  taken  to  be  a  sum  of  the 
Individual  fault  surface  rectangles,  each  being  weighted 
separately.  Unfortunately,  these  authors  used  a  rather 
simple  source  model  representation  in  that  it  »iad  dii  fer- 
ential  movement  only  on  a  horizontal  fault  parallel  to  the 
surface.  Furthermore,  they  gave  no  indication  as  to  how 
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they  arrived  at  their  final  model.  One  would  guess  that 
they  used  the  trial  and  error  method. 

Savage  and  Hastie  (1969)  and  Hastie  and  Savage  (1970) 
have  described  a  quasi-inversion  process  to  be  applied  to 
the  fitting  of  earthquake  static  displacement  data  using 
the  dislocation  models  for  an  arbitrarily  oriented  finite 
fault  surface  imbedded  in  a  homogenous  half-space.  In 
these  studies,  these  authors  swept  through  predetermined 
sets  of  sensitive  fault  parameters  —  fault  width,  dip 
angle,  depth,  and  slip  —  calculating  the  degree  of  fit  to 
all  the  data  for  each  model  tested.  he  model  which  best 
fit  the  data  in  a  least-squares  sense  was  termed  the  opti¬ 
mum  model.  These  calculations  seem  to  closely  coincide 
with  the  Monte  Carlo  techniques  used  to  find  acceptable 
models  of  the  radial  distributions  of  the  elastic  param¬ 
eters  within  the  earth  a3  described  by  Press  (1968,  1970, 
1972).  In  these  cases  a  reasonaole  fit  to  the  data  was 
obtained,  especially  in  the  case  for  the  Pairview,  Nevada 
earthquake.  Fitch  and  Scholtz  (1971)  later  extended  this 
work  to  some  degree.  However,  the  dislocation  model  used 
in  these  cases  was  highly  idealized  in  that  it  was 
restricted  to  the  Volterra  type  dislocation  in  which  the 
slip  was  constant  over  the  entire  dislocation  surface. 
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2 . 4  Construction  of  the  Pseudo-Somlgllana  Dislocation 

Static  Model  . 

Since  It  has  been  shown  that  the  fault  geometries  can 
be  more  complicated  than  just  plane  rectangular  surfaces, 
some  means  must  be  derived  to  allow  in  our  mathematical 
representation  of  the  faulting  process  for  these  complica¬ 
tions.  Complications  to  the  simple  models  can  occur  in  at 
leas  two  ways.  The  first  complication  is  that  we  wish  to 
be  able  to  allow  the  dislocation  to  take  on  arbitrary 
values  as  a  function  of  position  over  the  fault  surface. 
Secondly,  the  fault  surface  may  not  be  a  single  rectangular 
plane.  Both  of  these  complications  can  easily  be  repre¬ 
sented  approximately  by  discretizing  the  dislocation  sur¬ 
face.  That  is,  we  want  to  approximate  a  curved  fault 
surface  by  a  series  of  planar  surfaces  juxaposed  in  such  a 
manner  as  to  approximate  the  curvature  of  the  surface  to  be 
matched.  Curvature,  or  splaying,  could  be  thus  modeled  in 
any  direction.  An  example  of  matching  curvature  in  the 
horizontal  direction  could  be  envisioned  by  a  model  of  the 
San  Andreas  fault  which  includes  the  region  of  the  bend  in 
southern  California.  Here  a  series  of  plane  vertical  rec¬ 
tangular  surfaces  could  be  concatenated  horizontally  to 
match  the  observed  curvature.  Similarly,  a  dipping  thrust 
fault  in  which  the  dip  varies  with  depth  could  be  approxi¬ 
mated  by  a  series  of  rectangular  sheets  positioned  verti¬ 
cally  to  make  a  continuous  surface  in  which  the  dip  could 


change  discontinuously  between  fault  elements.  Examples  cf 
modeling  dipping  thrust  faults  in  this  manner  is  given  in 
later  chapters. 

With  this  same  scheme,  the  dislocation  could  be 
allowed  to  differ  on  each  of  the  surface  elements  which 
comprise  the  total  dislocation  surface.  Restriction  on 
the  variance  of  the  source  parameters  from  one  surface 
element  to  the  next  would  have  to  be  imposed  to  keep  the 
problem  physical. 

2.5  Linearization  of  the  Forward  Static  Problem. 

The  net  displacement  or  strain  field  at  the  surface, 
or  at  any  point  off  the  dislocation  surface  could  be  cal¬ 
culated  separately  for  each  of  the  individual  segments 
using  one  of  the  forward  problem  formulations  discussed 
earlier  in  this  chapter.  The  total  elastostatic  field  at 
a  particular  observation  point  would  be  a  simple  sum  of  the 
individual  contributions  from  each  of  the  comprising 
elements  . 

Wc  wish  to  pose  the  problem  in  such  a  way  as  to  be 
able  to  write  down  a  succint  relationship  between  the 
values  of  the  source  parameters  and  the  data  functionals 
which  we  compute  from  the  forward  problem  calculations. 
Suppose  that  we  calculate  the  values  of  the  elastostatic 
field  at  a  single  point  exterior  to  the  dislocation  sur¬ 
face  cf  our  chosen  fault  model  system  which  is  made  up  of 
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M  different  variously-oriented  dislocation  surface  ele¬ 
ments.  Consider  that  the  elastostatic  field  can  be 
described  by  N  field  variables,  preferably  those  for  which 
we  can  observe  in  the  field  following  an  earthquake. 
Suppose  that  there  are  L  source  parameters  which  can  be 
linearly  related  to  the  elastostatic  field  through  the 
forward  problem  formulations.  Then  this  relationship  is 
given  through  the  system  of  linear  equations 


L-M 

Y  A,.m,  =  d1  for  i  *  1,N.  (2.3) 

J-i  3  3 


In  these  sets  of  equations  d^^  are  the  calculated  elasto¬ 
static  field  functional  values,  m^  are  the  values  of  the 
linear  model  source  parameters,  and  the  coefficients  A^j 
are  the  elastic  media  response  of  a  particular  data  func¬ 
tional  due  to  a  particular  fault  surface  element  having  a 
unitary  source  strength  for  the  linear  parameters.  These 
coefficients  are  in  general  a  function  of  position.  If  we 
treat  the  components  of  m^  and  d^  as  elements  of  a  column 
and  row  vector  respectively  and  if  we  put  the  coefficients, 
j ,  in  standard  matrix  form  where  the  matrix  has  L*M 
columns  and  N  rows,  we  can  express  (2.3)  in  the  following 
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matrix  notation, 

Am  =  d  .  (2.4) 

The  model  components  of  m  are  contained  in  the  vector 
space  EL'M  and  the  data  functional  components  are  con¬ 
tained  within  the  vector  space  E^.  The  matrix,  A  ,  can  be 
considered  a  vector  operator  which  maps  E^  ^  into  E^ . 

We  have  been  careful  in  this  construction  to  limit 
ourselves  to  problems  where  the  source  parameters  in  the 
space  EL'M  can  be  linearly  related  to  the  calculated 
elastostatic  field  functionals  in  space  E^.  This  strictly 
linear  relation  is  valid  for  only  a  few  source  parameters 
in  special  instances.  If  the  forward  problem  is  to  be 
solved  by  the  analytic  closed  form  Green’s  function  solu¬ 
tions,  for  example  equation  (2.2),  then  we  have  to  impose 
the  Volterra  restriction 

Aui(P)  *  constant. 

With  this  restriction  we  can  write 

V«>  -  si*/^wVp’Q)vjdI  . 


and  the  problem  is  now  linear  with  respect  to  slip  in  the 


-29- 


ith  direction  on  the  individual  dislocation  surfaces.  In 
general,  the  solution  cannot  be  so  easily  linearized  with 
respect  to  other  parameters  which  characterize  the  dis¬ 
location  source-fault  length,  dip  angle,  depth,  position, 
etc..  An  examination  of  the  forward  equations  given  by 
Mansinha  and  Smylie  (1971)  is  convincing  with  this  respect. 
Fortunately,  by  numerically  evaluating  these  expressions, 
we  can  show  that  they  are  locally  linear.  The  extent  of 
the  locally  linear  domain  varies  from  source  parameters  to 
source  parameter  and  also  with  the  absolute  value  of  the 
source  parameter.  If  sufficient  care  can  be  paid  to  these 
details,  the  problem  can  be  approximately  linearized  for 
all  the  source  parameters  listed  above.  The  linearization 
can  be  accomplished  in  the  following  simple  way. 

The  degree  of  linearity  or  non-linearity  of  the 
forward  problem  functionals  for  the  various  source  param¬ 
eters  will  be  model  dependent,  that  is,  it  will  vary  from 
source  model  to  source  model.  If  we  wish  to  describe  the 
linear  domain  in  a  field  about  some  chosen  model,  n^,  we 
choose  some  other  source  model,  mg,  "near"  m1  such  that 
the  following  equation  can  be  written 

A  6m  =  fid  +  0[|m2  -m1|2].  (2.5) 
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The  following  definitions  have  been  applied: 

6m  «  m2  -  (2.6) 

6d  =  d(m2)  -  d(m^)  •  (2-7) 

d(m,)  is  the  elastostatic  field  for  a  particular  source 
model  m^.  The  problem  is  linearized  only  if  6m  is 
sufficiently  small  for  equation  (2.5)  to  hold.  The  con¬ 
ditions  for  linearity  discussed  here  are  equivalent  to 
requiring  the  forward  problem  functionals  to  be  Frechet 
differentiable  with  respect  to  the  source  parameters. 

If  we  calculate  the  forward  problem  for  a  source 
model  description  which  we  think  will  reasonably  approxi¬ 
mate  the  observed  static  field  functionals,  call  this 

model  m  ,  then  for  small  perturbations  about  this  model, 
s 

6n^,  an  approximate  linear  relationship  between  the  two 
vector  spaces  is  established.  Ihis  is  to  say  that  the 
coefficients  of  A^  are  linear.  We  note  here  that  in 
general,  the  coefficients  of  are  not  independent  of 
the  model  nig.  Indeed,  their  dependence  is  a  measure  of 
the  non-linearity  of  the  operator  coefficients  in  the 
region  of  the  model  space  being  sampled  by  the  test  model 
ni  .  The  perturbations,  dm  ,  must  remain  sir  all  in  the  sense 
that  they  are  approximately  linear  throughout  this  region. 
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In  our  matrix  notation  the  forward  problem  is  now  written 


A  6m  =  <5d  . 

S  o 


(2.8) 


2.6  Derivation  of  the  Stochastic  Inversion  Operators. 
Introduction .  This  section  addresses  the  problem  of  ob¬ 
taining  the  best  estimate  of  the  source  parameters  charac¬ 
terizing  a  fault  model  given  a  suite  of  observations  which 
can  be  linearly  related  to  the  faulting  process.  The 
problem  here  follows  closely  that  encountered  in  the 
studies  regarding  the  estimation  of  the  radial  distribu¬ 
tion  of  velocity  and  density  within  the  earth.  In  this 
area  of  research,  much  theoretical  progress  has  been  made 
in  the  last  six  years  in  the  treatment  of  inversion 
schemes  to  estimate  these  distributions.  Perhaps  the  most 
successful  and  certainly  the  most  elegant  of  these  schemes 
falls  in  the  general  category  of  stochastic  inversion 
theory.  This  theory,  which  will  be  applied  to  the  treat¬ 
ment  of  elastostatic  problems  in  this  thesis,  attempts  to 
give  the  best  estimate  of  a  discretized  approximation  to 
the  continuous  faulting  process  when  a  limited  amount  of 
data  is  obtainable.  As  pointed  out  by  Jordan  (1972),  the 
inverse  problem  when  posed  in  this  manner  usually  has  no 
unique  solution.  However,  the  solution  that  is  obtained 
is  unique  in  certain  respects,  as  will  be  discussed  later. 


-32- 


Furthermore,  the  stochastic  approach  allows  fo  the  in¬ 
clusion  of  inaccuracies  in  the  estimation  of  the  elasto- 
static  field  observations.  How  these  inaccuracies  affect 
our  model  estimations  will  be  fully  explored  in  the 
chapters  devoted  to  the  application  of  this  theory . 

The  fundamentals  of  the  theory  for  the  solution  of 
the  underconstrained  linear  inverse  problem  for  data  that 
contain  certain  amounts  of  "noise"  have  been  presented  by 
Backus  and  Gilbert  (1967,  1968,  1969,  1970).  Jordan  and 
Minster  (1971)  and  Jordan  (1972)  incorporated  portions  of 
the  Backus-Gilbert  theory  with  the  purely  stochastic 
theory  of  Franklin  (1970)  to  present  a  quite  complete 
approach  to  the  solution  of  this  type  of  problem.  The 
theory  as  applied  here  to  static  problems  is  essentially 
that  due  to  Jordan  (1972),  and  an  attempt  has  been  made 
to  follow  his  notation  throughout  this  thesis.  Sophisti¬ 
cated  discussions  as  to  the  validity  of  this  type  of  in¬ 
verse  and  the  mapping  functions  of  the  operators  are  given 
in  this  reference.  The  derivation  of  the  stochastic 
inversion  operators  below  are  given  only  in  the  context  as 
to  how  they  apply  to  the  elastostatic  problem.  In  the 
derivations,  for  reasons  of  simplicity  the  notation  used 
is  for  a  linear  problem.  If  applied  to  non-linear  problems 
that  have  been  linearized  in  the  procedure  discussed  above, 
the  difference  vectors  defined  in  (2.6)  and  (2.7)  are 
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merely  substituted  for  the  true  model  and  elastostatic 
field  vectors. 

Derivation.  Consider  the  problem  of  determining  some  M- 

M 

dimensional  vector  model,  m,  contained  in  the  space  E  , 
given  a  N-dimensional  elastostatic  field  vector,  d ,  in  the 
space  EN.  The  elastostatic  field  values  are  related  to 
by  the  system  of  linear  equations 

EAlJmJ  "  dl  i  •  1.  N  .  (2.9) 

J=1 

In  matrix  notation 

Ante  d  ,  (2.10) 

where  the  operator  A  solves  the  forward  problem  for  each  of 
the  N  elastostatic  field  values  contained  in  d  by  mapping 
EM  into  E^.  Thus  for  every  model  m  there  exists  some 
unique  determination  of  d  where 

d  -  d  (m)  .  (2.11) 

If  we  take  the  actual  field  observations  which  are  measured 
following  the  occurrence  of  an  earthquake  to  be  in  vector 
form,  do»  and  these  measurements  are  made  perfectly  with  no 


inaccuracies,  then 


-34- 


dQ  =  d  (m)  .  (2.12) 

Assuming,  of  course,  that  the  formulation  of  the  forward 
problem  will  exactly  determine  the  elastostatic  field 
values.  However,  if  there  are  any  inaccuracies  in  the 
observed  field  values,  then  these  observed  values,  dQ,  can 
be  written  as  a  combination  of  the  projected  field  values 
plus  some  measure  of  the  uncertainty  in  these  observations, 

dQ  =  d(m)  +  n  ,  (2.13) 

or  by  substitution  from  (2.10) 

A  m ■  dQ  +  n  .  (2 . 14) 

Here  n  is  a  vector  containing  the  components  of  the  "noise" 

in  the  observed  field  data.  We  assume  that  this  noise  is 

randomly  distributed  in  a  Gaussian  fashion  and  that  any 

bias  to  the  data  is  removed  before  the  noise  is  estimated. 

Each  component,  dn  ,  is  assumed  to  be  the  mean  of  a 

1  2 

Gaussian  random  variable  with  variance,  a^.  We  can  define 
a  diagonal  variance  operator,  Cnn,  to  be 
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(2.15) 


where  is  the  estimated  variance  of  the  ith  data  value. 
In  assuming  this  diagonal  form,  we  are  implicitly  assuming 
that  there  is  no  co-variance  between  data. 

Since  the  only  information  that  we  have  about  m  is 
contained  in  (2.14) ,  we  know  nothing  about  the  components 
of  m  which  lie  outside  the  space  RCEM  which  is  spanned  by 
the  vectors  {a.,:i*l,N}.  It  is  reasonable  to  require  that 

a- 

our  estimate  of  m,  call  it  m,  lie  totally  within  the  sub¬ 
space  R;  then  we  can  assign  a  non-zero  value  to  only  those 
components  for  which  we  have  information.  Under  this 
restriction,  we  can  write 


m  = 


A*  b 


(2.16) 


for  some  vector  b  contained  in  the  vector  data  space  E  . 

A* 

In  this  last  equation,  we  are  using  the  notation  A  to 
represent  the  transpose  of  the  matrix  A.  This  convention 
will  be  used  throughout  this  thesis.  To  select  an  optimal 
b,  call  this  b,  we  wish  to  minimize  a  suitable  quadratic 
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measure  of  the  errors  involved  in  this  estimation.  We 
choose  to  minimize  some  weighted  sum  of  two  measures  of 
the  errors  involved  in  this  problem.  This  weighted  sum  can 
be  parameterized  by  a  trade-off  curve  between  these  two 
errors,  with  the  position  along  this  curve  used  as  the 
parametric  factor.  Specifically,  we  want  to  minimize 

e2(0,b)  =  e2(b)  cos(6)  +  e2(b)  sin(0)  ,  (2.17) 


where 


e2(b)  -  |  |m  -  A*b|  |  2  f  (2.18) 

and 

r2(b)  -  i  Cnn  b  ,  (2.19) 

p 

The  first  measure  of  error,  e1(b),  is  the  square  of  the 
Euclidian  norm,  defined  by 

?  Mo 

II  X  ||2  -  £  x2  , 

i-1 

of  the  difference  between  our  estimate  of  the  model,  m, 
and  the  actual  vector  we  are  estimating.  This  quantity 
decreases  as  we  more  closely  approximate  m.  The  second 
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measure  of  error  associated  with  b  arises  from  uncertainty 
in  the  components  of  dQ.  This  quantity  decreases  as  our 
estimate  becomes  more  reliable.  The  parameterization 
angle,  6,  is  allowed  to  vary  on  Che  interval  [0,tt/2],  so 
that  at  6  =  ir/2,  e^Cb)  is  minimized,  indicating  maximum 
reliability  of  the  model.  At  0  =  0,  e^(b)  is  minimized, 
indicating  maximum  accuracy  in  the  estimation  of  the  model. 

We  note  here  that  these  two  errors  are  measured  with 
two  different  norms,  each  in  the  model  space.  We  must 
establish  some  common  norm  on  each  of  these  errors  so  that 
the  parameterization  of  the  sum  of  these  errors  can  be 
accomplished.  This  normalization  is  performed  through  the 
introduction  of  a  correlation  operator, W.  This  correla¬ 
tion  operator  can  be  thought  of  simply  in  terms  of  a 
weighting  function  for  the  various  model  components .  The 
norm  of  this  operator  is  fixed  so  that  at  the  critical 
point  on  the  trade-off  curve  between  the  two  types  of 
errors,  at  6  «  ir/i* ,  the  absolute  value  of  the  two  errors 
are  equal . 

For  the  present,  we  assume  that  the  correlation 
operator,  W,  is  the  idemfactor,  I  ,  so  that  this  effect  can 
be  ignored  in  our  minimization  calculations.  The  results 
of  this  minimization  then  will  be  generalized  to  include 
an  arbitrary  correlation  operator. 
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In  order  to  minimize  e2(8,b)  given  in  equation  (2.17) > 
we  take  6b  to  be  a  small  arbitrary  perturbation  of  b.  To 
first  order  in  6b  we  can  write 

6e2(6,b)  =  e2(0,b  +  6b)  -  e2(0,b)  . 

Performing  this  first  order  perturbation  on  equation  (2.17) 
we  find  that 

6e2(0 ,b)  =  2[bAA*-  m/f ]6b  cos(0)  +  2  bCnn6b  sin(0)  . 

In  order  to  minimize  e2(0,b),  we  set  6e(0,b)  =  0.  When 
this  is  done,  we  see  that  fib  truly  is  an  arbitrary  pertur¬ 
bation,  and  e2(0 ,b)  will  be  stationary  if  and  only  if 

(A  A  +  tan(  6)  Cnn)b  -  Am  .  (2.20) 

It  can  be  shown  (Jordan,  1972)  that  this  stationary  point 
is  a  unique  minimum,  and  the  vector,  b  ,  which  satisfies 
this  condition  is  our  optimum  vector,  b. 

p 

If  C  is  non-singular,  that  is,  each  o,  ¥  0,  and 
nn  x 

6  >  0,  then  the  matrix  (AA  +  tan(0)  Cnn)  is  non-singular 
and 

b«  (AA*  +  tan(6)  Cnn)_1Am  • 


(2.21) 
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In  this  last  equation,  m  is  unknown  but  by  substituting 
from  equation  (2.10)  we  get 


b=  (  AA*  +  tan(e)  Cnn)_1  d0  • 


(2.22) 


Substituting  this  optimal  value  of  b  into  equation  (2.16) 
we  see  that  the  optimal  estimate  of  the  model  for  a  fixed 
value  of  9  will  be  given  by 

m  -  A*  (AA*  +  tan(e)  C^)'1  d0  .  (2-23) 


In  the  above  results  all  components  of  m  are  equally 
weighted  with  the  identity  operator.  A  more  general 
weighting  can  be  introduced  by  considering  a  set  {Wjjj*l,M} 
of  non-zero  positive  weights  for  the  model  components. 

Let  us  define  this  weighting,  or  correlation  matrix,  in  the 
following  manner, 
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This  leads  us  to  define  a  normalized  elastic  media  response 
operator 

A'  =  AW  .  (2.24) 

With  this  normalized  definition,  equation  (2.10)  is  now 
written 

A W-1  m  ■  d0  .  (2.25) 

Following  the  same  procedure  as  before,  we  require 

m  *  A  b  ,  (2.26) 

and  minimize 

t 

ey( 6 ,b)  *  | | m -  A  *b| |y  cos(6)  +  lfCnnb  sin(0)  , 
where  ||’||w  is  the  weighted  norm  defined  by 

mu*  ■  z’w'lz  ■  • 

This  weighted  norm,  of  course,  reduces  to  the  Euclidian 

norm  if  -  1  for  all  i  -  1,M.  The  minimization  of 

e2(0,b)  with  respect  to  a  variation  of  b  procedes  as  before. 
W 

The  results  are 

m«  WA(AWA*  +  tan(e)  C^)"1  dQ 


(2.27) 


-41- 


Now  since  there  are  uncertainties  in  the  observed 
elastostatic  field,  the  best  estimate  of  the  model,  m,  is 
some  filtered  average  of  the  true  model,  m,  given  by 

m  =  R  m  (2.28) 

This  averaging  operator,  which  contains  the  response  ker¬ 
nels  for  the  elements  of  mean  easily  be  found  by  substi¬ 
tuting  for  dQ  in  equation  (2.25).  Performing  this  substi¬ 
tution  in  equation  (2.27)  we  obtain 

m  =  WA'(AWA*  +  tan(e)  Cnn)‘]-A  m  , 
or  by  inspection  from  equation  (2.28) 

R  -  WA(AWA*  +  tan(0)  Cnn)_1A  .  (2-29) 

Individual  rows  of  this  operator  contain  the  averaging 
of  the  estimated  values  of  the  Individual  model  components 
with  respect  to  the  other  model  components.  This  averaging 
is  taking  place  in  a  sense  that  the  estimation  of  the  ith 
model  component  is  actually  the  true  value  of  this  component 
’’convolved’'  in  the  model  space  with  the  function  defined  on 
the  model  space  by  the  components  of  the  ith  row  of  the 
averaging  operator.  If  a  particular  model  component  is 
perfectly  determined,  say  the  ith  value,  that  is,  its  value 
is  perfectly  resolvable,  then  *  1  and  all  other  -  0. 
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In  the  limit  of  infinite  resolution  on  all  model  com¬ 
ponents,  that  is,  either  0  “  0  or  Cnn  =  0,  the  averaging 
operator  approaches  the  idemfactor,  I  . 

By  similar  substitutions,  we  can  express  equation 

(2.19)  as 


This  operator  is  termed  the  variance  operator.  Examining 

equation  (2.30)  we  see  that  the  bilinear  product  of  this 

operator  and  the  model  components  is  a  measure  of  the  error 

induced  from  the  data  space,  through  the  variance  matrix 

C  .  into  the  model  space .  Since  we  are  assuming  that  the 
nn 

errors  exhibited  in  Cm  are  normally  distributed,  we  can 
determine  the  following  about  the  errors  induced  from  the 
data  space  due  to  inaccuracies  in  the  description  of  the 
elastostatic  field  into  inaccuracies  in  the  estimated 
source  model  parameters.  Use  of  this  operator  does  not 
tell  us  the  absolute  Inaccuracies  of  our  estimated  model 
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per  se ;  instead,  it  can  only  tell  us  whether  or  not  a 

certain  perturbation  in  the  model  is  resolvable  to  a  cer¬ 
tain  degree  by  the  data.  So  in  practice,  we  have  to 
prescribe  a  perturbation  vector  on  our  source  model  and  test 
to  see  if  the  data  can  "see"  thi3  perturbation.  This 
ability  to  distinguish  model  perturbations  by  the  observed 
data  will  depend  directly  on  the  accuracy  of  the  data.  The 
more  accurate  the  data,  the  smaller  a  model  perturbation 
these  data  will  be  able  to  detect.  Since  we  are  now  map¬ 
ping  errors  in  the  opposite  direction  as  that  defined  in 
equation  (2.30),  clearly  the  inverse  of  this  operator  is 
the  projection  that  we  desire.  Since  the  errors  are 
induced  in  directions  along  the  eigenvectors  of  V(6),  then 
we  choose  to  take  the  inverse  of  this  operator  as  the 
generalized  inverse  given  by 

v+(e>  •  (2'32) 

i-1  Xi 

Here  we  are  assuming  that  V(0)  hap  a  total  of  J  non-zero 
eigenvalues  ( A ^ ,  i«l,J)  with  the  associated  eigenvectors  u1> 

ft 

The  notation  u1«u1  indicates  an  outer-product  expansion 

'  *  + 

between  the  two  vectors  and  u^.  Since  V  (0)  is  a 
generalized  inverse  of  V(0),  then  the  inner  product  of 
V(0)  and  V+(0)  are  not  necessarily  the  Identity  operator 
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but  rather  some  projection  operator,  Pv  that  is  both 
idempotent  (P  Py  =  Py)  and  symmetric  (Py  =  Py)  . 

In  particular,  some  vector  perturbation  in  the  model 
space ,  q  ,  is  resolvable  to  within  a  certain  confidence 
limit,  to  which  we  can  assign  some  confidence  coefficient 
k(c),  if  the  following  inequality  holds. 

q*  V*  Q  >  k2(c)  •  (2.33) 


For  example,  for  the  95^  confidence  limit,  k(c)  can  be 
found  in  any  good  statistics  reference  to  be  1.96. 

A  two-dimensional  geometrical  argument  will  illus¬ 
trate  the  use  of  equation  (2.33)*  Assume  that  the  errors 

induced  from  the  data  space  onto  the  model  space  by  the 

2 

variance  operator  (eigenvalues  of  this  operator)  are  om 

p  ^ 

and  a  .  (This  variance  should  not  be  confused  with  the 
in^i 

data  variance  defined  in  equation  (2.15)).  These  errors 

A  A 

lie  along  the  eigenvector  directions,  x^  and  x2  respect¬ 
ively.  Now  if  a  vector  x  has  components  along  these 
directions  then  the  equation 


x*v\ 


can  be  written  out 


m. 


_2_ 

’L 
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This  is  Just  the  equation  of  an  ellipse  whose  semi-major 

axes  are  ko  and  kom  .  This  ellipse,  or  hyper-ellipsoid 
m  ^  iti^ 

when  this  argument  is  extended  to  higher  dimensions,  is 
called  the  confidence  ellipse.  The  enclosure  of  this 
ellipse  represents  the  area  of  unresolvable  model  pertur¬ 
bations,  and  the  area  exterior  to  the  ellipse  represents  a 
model  perturbation  which  is  large  enough  to  be  resolvable 
by  the  data  at  a  certain  confidence  limit  associated  with 
the  axis  parameter  k.  By  making  k  larger,  we  are  increas¬ 
ing  the  confidence  limit  and  increasing  the  size  of  the 
confidence  ellipse  thus  requiring  larger  model  perturba¬ 
tions  before  they  can  be  detected  by  the  data  at  that  con¬ 
fidence  limit.  In  order  to  check  the  resolvability  of  a 
given  model  perturbation, we  choose  our  value  of  k  (say 
1.96)  and  merely  test  to  see  if  this  vector  protrudes  the 
confidence  ellipse.  We  note  here  that  this  resolvability 
criterion  depends  only  on  relative  perturbations  to  the 
source  model  parameters  and  not  on  the  absolute  configura¬ 
tion  of  the  final  or  optimum  model  that  we  obtain  from  the 
inversion  process.  Thus  we  have  to  propose  a  hypothetical 
perturbation, or  a  series  of  perturbations.  Judiciously 
chosen  to  explain  or  disclaim  certain  features  of  our 
model,  and  expose  them  to  this  testing  procedure.  Only  on 
this  basis  can  we  determine  the  limits  the  model  source 


-116- 


parameters  can  take  and  still  fit  the  observed  ela.3to- 
static  field.  The  power  of  this  operator  becomes  apparent 
when  applied  to  actual  problems  as  we  shall  see  in  later 
chapters . 

2 . 7  Discussion . 

In  this  chapter  we  have  discussed  the  development  of 
methods  of  obtaining  an  accurate  representation  of  the 
forward  elastostatic  problem  for  a  given  description  of 
the  faulting  process.  We  have  reviewed  the  early  uses  of 
these  forward  formulations  in  attempting  to  deduce  source 
parameters  which  can  characterize  a  given  event.  A  method 
was  suggested  by  which  a  more  complicated  and  arbitrary 
static  dislocation  function  could  be  approximated  with  the 
formulations  derived  from  simple  dislocation  sources.  It 
was  found  that  by  making  possible  a  more  complex  static 
source  description  some  means  must  be  used  to  systemati¬ 
cally  relate  the  observed  elastostatic  phenomena  to  the 
media  response  from  f  *  various  source  parameters.  The 
stochastic  inversion  scheme  provided  an  ideal  means  to 
give  the  best  estimates  to  the  solution  for  the  usually 
underdetermined  static  problem.  By  use  of  this  inversion 
scheme,  we  can  benefit  from  the  use  and  knowledge  of  the 
various  operators  which  fall  out  of  the  derivations.  These 
operators  deal  with  the  errors  in  both  the  observations  and 
those  in  our  solutions.  Quantitative  appraisals  of  the 
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decency  of  a  given  solution  to  a  particular  static  problem 
become  available  through  the  use  of  these  operators. 

For  the  special  case  of  0  =  0,  equation  (2.23)  is 
commonly  known  as  the  generalized  inverse.  For  this  case. 
Noble  (1969,  p.  143)  has  shown  through  the  method  of 
Lagrange  multipliers  that  the  generalized  inverse  also 
minimizes  the  norm  of  m.  We  can  think  of  this  as  physi¬ 
cally  giving  the  longest  wavelength,  or  smoothest  model 
solution,  for  a  given  set  of  data.  In  elastostatic  prob¬ 
lems,  this  property  is  especially  valuable,  since  we  would 
expect  the  displacement  on  a  fault  surface  to  locally  vary 
in  some  fairly  smooth  fashion. 

By  combining  all  of  the  formalisms  discussed  in  this 
chapter,  we  should  be  able  to  take  a  formidable  advance  in 
our  understanding  of  the  static  processes  which  accompany 
earthquakes.  The  theory  discussed  here  will  be  applied  to 
data  from  actual  earthquakes  in  the  following  chapters. 
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Chapter  3 

A  Static  Dislocation  Model  of  the 
1964  Alaska  Earthquake 


3.1  Introduction . 

The  Alaska  earthquake  of  28  March  1964  which  was 
centered  near  Prince  William  Sound  was  probably  the  largest 
seismic  event  in  North  America  this  century.  The  magni¬ 
tude  of  this  event  has  been  estimated  to  be  between 

M  *  8.3  to  M  =  8.6.  With  the  possible  exception  of  the 
s  s 

1971  San  Fernando,  California  earthquake,  this  earthquake 
has  been  the  most  intensely  studied  occurrence  in  the 
history  of  geophysics.  The  regional  deformation  accompany¬ 
ing  this  event  involved  changes  In  land  level  of  unprece- 

2 

dented  areal  extent,  encompassing  some  200,000  km  .  The 
residual  vertical  displacements  produced  were  measurable 
geodetically  along  a  400  km  profile  approximately  perpen¬ 
dicular  to  the  Gulf  of  Alaska  and  approximately  800  km 
adjacent  and  parallel  to  the  coastline.  Yet  despite  the 
importance  that  this  earthquake  had  on  the  tectonic 
character  of  the  affected  region  and  the  importance  of  the 
contributions  that  the  data  from  this  event  provided  toward 
an  increased  scientific  understanding  of  the  origin  of 
earthquakes,  considerable  controversy  still  surrounds  the 
exact  source  mechanism.  It  is  hoped  that  the  results  from 
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this  chapter  will  help  allay  some  of  this  controversy. 

3 . 2  Fault  Representation. 

Since  the  first  studies  of  the  196M  Alaska  earthquake, 
the  main  focal  mechanism  and  the  accompanying  sense  of 
motion  have  remained  somewhat  of  a.  controversy  primarily 
because  of  the  ambiguity  of  the  fault  plane  solutions 
based  on  P-wave  first  arrival  data.  The  two  contesting 
mechanisms  are  one  having  the  geometry  of  a  nearly  verti¬ 
cal  reverse  fault,  and  the  other  a  low  angle  thrust  fault. 
Figure  3.1  shows  a  profile  extending  from  the  southeast  to 
the  northwest  approximately  bisecting  the  elongated  area 
of  deformation.  This  cross  section  corresponds  to  profile 
BB«  shown  in  Figure  3-2.  In  Figure  3-1  we  have  diagrammat- 
ically  represented  the  two  possible  fault  plane  mechanisms 
and  their  relation  to  the  hypocenter,  shown  at  the  inter¬ 
section  of  these  two  planes.  The  representative  geometry 
that  we  choose  to  explain  in  detail  the  static  fields 
which  accompanied  this  earthquake  must  be  in  reasonable 
compatibility  with  the  geometry  necessary  to  explain  the 
following  observed  or  calculated  entities: 

1)  epicenter  location  and  hypocentral  depth 

2)  P-wave  first  motion  polarities  and  S-wave 
polarizations 

3)  aftershock  distribution 

4)  radiation  patterns  of  long  period  Love  and 


Figure  3.1.  Schematic  diagram  of  the  two  possible  nodal 
planes  and  the  relative  dislocation  on  each.  The  hypocenter 
of  the  main  shock  Is  located  at  the  intersection  of  the 
two  planes. 


Figure  3.2.  Regional  deformation  that  accompanied  the 
March  28,  1964,  Alaska  earthquake.  Cross  section  used  in 
this  study  is  labeled  BB'.  The  Patton  Bay  fault  is 
indicated  as  the  axis  of  uplift.  Figure  is  from  Plafker 
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Ray  leigh  waves 

5)  geological  reasoning  for  faulting  —  island  arc 
implications 

6)  near  field  displacements. 

We  will  briefly  review  the  geophysical  literature  for 
supportive  arguments  to  favor  one  or  the  other  of  the  pro¬ 
posed  mechanisms.  We  will  then  adopt  a  model  which  we 
think  will  best  fit  all  of  these  criterior. 

The  hypocentral  depth  for  the  main  event  was  first 
given  to  be  about  20  km,  and  in  later  calculations  with 
the  inclusion  of  more  data  the  depth  was  restricted  to  33 
km.  (This  restricted  depth  is  the  standard  depth  assigned 
a  shallow  event  when  the  depth  determination  algorithm 
does  not  converge,  or  else  converges  to  a  negative  depth.) 
No  depth  sensitive  phases,  such  as  pP  or  sP  could  be 
positively  identified  on  records  of  the  main  shock.  A 
reasonable  assumption  would  be  to  place  the  depth  as  lying 
between  20  km  and  50  km.  The  hypocenter  certainly  was  not 
deep  as  evidenced  by  the  large  amplitude  surface  waves 
generated  by  this  earthquake.  The  epicenter  of  the  main 
shock  was  located  by  Sherburne  et_  al .  (1968)  and  von  Hake 
and  Cloud  (1966)  to  be  near  the  north  shore  of  the  Prince 
William  Sound  on  the  small  peninsula  separating  College 
Fiord  and  Unakwik  Inlet.  The  coordinates  of  the  epicenter 
are  given  as  61.0*4°  +  0.05°  north  latitude  and  1*47.73° 


+  0.07°  tfest  longitude. 

The  focal  mechanism  for  the  first  motion  of  this 
earthquake  has  been  studied  by  a  number  of  authors 
( Algermissen ,  1964,  1965,  1966;  Harding  and  Algermissen, 
1969;  Berg,  1965;  Stauder  and  Bollinger,  1966).  These 
studies  show  only  one  fairly  well-defined  nodal  plane. 

There  is  some  slight  ambiguity  in  the  exact  orientation  of 
this  plane  due  to  non-impulsive ,  or  emergent  P-wave  first 
arrivals  at  a  number  of  key  stations,  but  this  is  a  second 
order  effect.  The  preferred  orientation  of  this  nodal 
plane  is  given  to  be  strike  N  62°  E,  dip  82°  S.  The  defi¬ 
nition  of  the  second  nodal  plane  is  limited  because  of  the 
almost  total  lack  of  geographica:  control  in  the  station 
locations.  Berg  (1965)  attempted  with  limited  success  to 
determine  the  orientation  of  this  second  nodal  plane  by 
observing  a  dilatation  at  one  station,  Yellowknife,  Canada. 
The  location  of  this  station  is  critical  in  defining  this 
second  nodal  plane.  The  orientation  of  this  plane  has  been 
estimated  to  give  a  dip  of  26°  to  the  northeast.  This 
unfortunate  distribution  of  stations  to  the  north  of  the 
epicenter  precludes  the  identification  of  the  nodal  plane 
that  would  be  present  due  to  a  low  angle  thrust,  although 
the  plane  has  been  restricted  by  the  data  presented  by 
Stauder  and  Bollinger  (1966).  These  authors  conclude  that 
the  second  nodal  plane  can  have  a  dip  varying  from  less 
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than  25°  to  the  northeast  through  5°  to  the  northwest  to 
less  than  60°  .  '  the  southeast.  S-wave  polarization  studies 
suffer  from  the  sa^e  restriction  in  the  station  distribu¬ 
tion  respect.  The  results  from  the  S-wave  polarization 
angle  study  by  Harding  and  Algermissen  (1966)  indicate  that 
for  a  double  couple  type  source  on  a  nearly  vertical  fault 
the  required  motion  to  fit  the  observed  S-wave  polarities 
would  be  predominantly  strike-slip. 

One  suggestion  that  must  be  kept  in  mind  when  trying 
to  interpret  the  orientation  of  the  nodal  planes  from  first 
motion  data  is  that  presented  by  Wyss  and  Brune  (1967)  . 

These  authors  suggested  that  the  faulting  which  occurred 
over  the  entire  segment  involved  a  complex  multiple  rupture 
mechanism.  If  this  mechanism  is  in  fact  the  way  the  fault¬ 
ing  took  place,  then  the  Initial  motion  at  the  hypocenter 
can  have  little,  if  any,  bearing  on  how  the  faulting  pro¬ 
ceeded  as  a  whole . 

One  clue  as  to  the  possibility  of  deciding  which  type 
faulting  took  place  is  given  by  examining  the  spatial  dis¬ 
tribution  of  aftershocks.  Algermissen  e_t  al .  (1972)  present 
just  such  data  for  over  2,000  locatable  aftershocks. 

Special  attention  was  given  to  a  sub-set  of  this  aftershock 
location  data  which  were  well  located  and  contained  posi¬ 
tively  Identifiable  depth  phases.  These  events  showed  that, 
especially  in  the  vicinity  of  Prince  William  Sound,  the 


-55- 


aftershocks  were  shallow.  In  fact,  approximately  62$ 
were  located  at  depths  less  than  20  km  with  only  1$  of  the 
events  located  at  depths  greater  than  40  km.  This  depth 
distribution  of  aftershocks  suggests  that  most  or  all  of 
the  faulting  was  confined  within  the  crust  and  perhaps  the 
top  of  the  upper  mantle  along  the  continental  margin. 

These  authors  depict  the  foci  of  the  aftershocks  located 
in  this  area  under  consideration  as  defining  a  plane  which 
dips  at  a  shallow  angle  (4°-6°)  under  the  continental  block. 
Focal  mechanism  studies  of  the  aftershock  by  Stauder  and 
Bollinger  (1966)  delineate  a  fault  plane  some  600  km  in 
length  and  at  least  200  km  in  width  having  an  average  dip 
of  about  10°,  while  the  main  shock  had  a  depth  of  focus  of 
between  20  and  50  km  and  had  a  body  wave  nodal  plane  solu¬ 
tion  dipping  between  10°  and  15° . 

The  outer  limits  of  the  aftershock  region  appear  to  be 
very  well  defined  and  the  region  is  not  confined  along  the 
surface  trace  of  the  postulated  steep-fault  model.  The 
aftershocks  lie  mainly  in  a  bread  belt  roughly  paralleling 
the  continental  margin  mostly  falling  in  the  area  of  mapped 
or  inferred  major  uplifting.  The  aftershock  zone  is  not 
even  approximately  centered  on  the  epicenter  of  the  main 
shock . 

Surface  wave  studies  of  this  earthquake  have  been 
limited  to  long  period  multiple  Love  and  Rayleigh  waves  due 
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to  the  tangled  complexity  of  the  large  amplitude  records 
at  the  WWSSN  stations.  With  only  few  exceptions,  che  first 
multiples  to  be  fully  recovered  have  been  the  R4  and  G4 
wave  trains .  These  signals  have  been  analyzed  In  two 
different,  but  hopefully  equivalent,  ways.  Toksoz  et  al . 
(1965)  and  Ben  Menahem  et  al.  (1972)  used  the  spectral 
phase  and  amplitude  equilization  method  while  Kanamori 
(1970)  used  a  time-domain  at.  ■  ysis.  For  a  simple  point 
double  couple  source,  the  radiation  patterns  for  surface 
waves  for  the  two  contesting  fault  orientations  are  approx¬ 
imately  equivalent.  However,  if  the  source  has  some 
finiteness  as  exhibited  by  propagating  in  a  given  direc¬ 
tion  then  assymmetries  in  the  Love  and  Rayleigh  wave  ra 
diation  patterns  are  introduced.  As  pointed  out  by  Savage 
and  Hastie  (1966,  p.  4899-4900),  the  assymmetries  between 
Love  and  Rayleigh  wave  radiation  patterns  will  be  different 
only  if  the  rupture  propagation  is  not  along  the  null 
axis.  If  rupture  does  take  place  in  a  direction  away  from 
this  axis  then  there  is  a  possibility  of  distinguishing 
uniquely  the  two  fault  orientations.  Because  of  differ¬ 
ences  in  azimuthal  cov  *age,  Ben  Menahem  et  al.  did  not 
detect  any  assymmetries  in  his  radiation  patterns  while 
Kanamori  did.  Kanamori  interprets  this  assymmetry  in 
terms  of  a  measured  component  of  rupture  propagation 
normal  to  the  strike  of  the  ^ault.  His  solutions  favor 
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the  low-angle  thrust  mechanism  and  his  model  is  compatible 
with  the  long  period  surface  waves  radiation  patterns  for 
a  fault  dipping  at  about  20°. 

Plafker  (1965)  uses  his  interpretation  of  a  vast 
quantity  of  field  observations  in  the  area  of  deformation 
to  argue  rather  forcefully  for  the  low  angle  thrust  mechan¬ 
ism.  These  arguments  will  not  be  repeated  here  but  are 
based  mainly  on  the  large  displacements  in  relation  to  the 
focal  mechanism  studies  and  the  spatial  distribution  of 
aftershock  seismicity.  Plafker  (1972)  extends  much  the 
same  arguments  for  a  low  angle  thrust  fault -in  the  context 
of  being  consistent  with  the  mechanism  expected  for  island 
arc  tectonics  (Isacks  et  al. ,  1968;  Stauder,  1968).  He 
concludes  that  the  earthquake  occurred  as  shear  failure  on 
a  fairly  complex  major  low  angle  thrust  fault,  or  mega¬ 
thrust,  that  dips  from  the  vicinity  of  the  offshore  trench 
to  beneath  the  continental  margin.  The  overthrusting  is 
interpreted  in  terms  of  elastic  rebound  resulting  from  the 
progressive  unde rthrus ting  of  the  oceanic  crust  and  mantle 
beneath  the  continental  margin  prior  to  1964.  This  mech¬ 
anism  is  consistent  with  Benioff's  (1954)  theory  for 
oceanic  trenches  and  associated  mountain  ranges. 

On  the  basis  of  modeling  the  observed  vertical  dis¬ 
placements,  Press  and  Jackson  (1965)  and  Press  (1965) 

\ 

attempted  to  demonstrate  that  the  observed  uplift  and 
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subsidence  could  be  accounted  for  by  about  10  m  of  con¬ 
stant  dip-slip  motion  on  a  vertical  plane  extending  from 
a  depth  of  about  15  km  down  to  a  depth  of  150  km  or  more. 
These  authors  did  not  include  in  their  data  set  all  verti¬ 
cal  displacement  points  available.  Savage  and  Hastie 
(1966)  and  Hastie  and  Savage  (1970)  got  better  results 
trying  to  fit  the  same  data  with  a  low  angle  thrust  fault 
with  about  10  m  of  constant  displacement  over  the  entire 
surface.  Savage  and  Hastie  showed  that  the  vertical  re¬ 
verse  fault  model  geometry  placed  the  zone  of  maximum  sub¬ 
sidence  too  close  to  the  zone  of  maximum  uplift,  whereas 
for  the  low  angle  thrust  geometry,  this  observed  lack  of 
symmetry  in  the  vertical  displacements  is  approximately 
satisfied.  Stauder  and  Bollinger  (1966)  accomplished  a 
more  realistic  modeling  of  the  displacements  on  a  hori¬ 
zontal  thrust  fault  on  which  differential  movement  on  the 
fault  surface  was  allowed.  These  authors  tried  to  include 
the  effects  of  local  or  subsidiary  faulting  on  Montague 
Island  (see  Figure  3*2).  The  lccal  faulting  shows  a 
dominance  of  vertical  slip  and  has  been  described  by 
Plafker  (1965)  and  Grantz  et  al^  (1964a, b).  Stauder  and 
Bollinger  (1966)  model  this  secondary  fault  as  a  constant 
dip-slip  dislocation  on  a  vertical  surface  directly 
beneath  Montague  Island. 

Additionally,  the  low  angle  geometry  is  preferable  in 


describing  the  behavior  of  the  observed  extensive  hori¬ 
zontal  surface  deformation  as  reported  by  Parkin  ( 1 96 6 ) . 

The  sense  of  this  deformation  is  mainly  consistent  with  the 
seaward  overthrusting  of  the  continental  block.  This 
direction  of  motion  is  especially  predominant  in  the  area 
between  the  Kenai  Mountains  and  the  offshore  islands. 
However,  we  see  from  Figure  3.1  that  we  would  intuitively 
expect  the  horizontal  displacements  to  be  in  the  opposite 
direction  if  the  steeply  dipping  reverse  faulting  mechanism 
were  adopted.  Thus  we  have  decided  to  adopt  the  low  angle 
thrust  geometry  for  our  fault  model  in  explaining  the 
surface  displacement  data  because  it  seems  most  consistent 
with  the  seismic,  geodetic  and  geologic  observations  per¬ 
taining  to  this  earthquake. 

In  each  of  these  attempts  in  modeling  the  vertical 
displacements  the  formulation  of  a  dislocation  in  a  uni¬ 
form  elastic  half-space  was  used  (Green's  functions 
solutions).  Since  this  is  a  region  where  there  is  a  large 
contrast  in  the  juxaposed  crustal  types  —  oceanic  crust 
underthrusting  continental  crust  —  this  uniform  elastic 
half-space  approximation  may  not  be  appropriate.  This 
approximation  will  be  investigated  later  in  this  chapter. 
All  of  the  above  models  are  able  to  fit  only  the  gross 
features  of  the  zero-frequency  data  of  this  earthquake, 
not  Just  because  the  earth's  crust  is  not  a  uniform  elastic 
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half-space  and  the  slip  varies  continuously  along  the 
fault  plane,  hut  also  because  the  estimates  of  the  fault 
offsets  were  not  related  to  the  observations  in  a  system¬ 
atic  fashion. 

For  this  earthquake  we  will  model  the  tectonic  en¬ 
vironment  with  a  laterally  heterogenous  geologic  model. 

The  finite-element  formulation  will  be  used  to  compute  the 
static  response  of  a  structural  model  of  the  crust  to  a 
unit  offset  imposed  on  a  series  of  nodal  segments  repre¬ 
senting  the  fault,  and  the  inversion  technique  will  be 
used  to  invert  any  free-surface  statical  observations  to 
obtain  the  proper  linear  combination  of  these  offsets  which 
will  result  in  a  computed  movement  of  the  surface  which 
fits  the  observed  data  to  some  chosen  degree  of  accuracy. 
Since  the  finite-element  formulation  used  in  this  chapter 
is  limited  to  solving  problems  involving  plane  strain 
elasticity,  any  displacement  profile  that  is  to  be  modeled 
correctly  must  be  approximately  free  of  fault  end  effects 
and  movement  due  to  strike  slip  motion.  The  effect  of 
assuming  an  infinite  length  fault  will  be  discussed  in  a 
later  portion  of  this  chapter. 

Tne  structural  model  chosen  for  this  svudy  is  given 
in  Figure  3.3.  The  geometry  is  based  upon  that  suggested 
by  Plafker  (1972)  and  Stanley  (1966)  as  being  the  mo  t 
consistent  with  the  regional  tectonic  setting  of  the 
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Patton  Bay  fault  are  indicated  by  heavy  lines.  The  epicenter  as  projected 
onto  this  profile  is  given  by  the  "star”.  The  arrows  on  the  lateral  bound 
aries  indicate  the  imposed  displacement  boundary  condition  used  to  model 
the  average  regional  prestress  field  expected  from  tectonic  plate  motion. 
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earthquake,  seismic  refraction  studies,  and  the  earthquake 
distribution  of  the  area  (Tobin  and  Sykes,  1966).  The 
region  is  modeled  by  four  geologic  provinces  ,  and  the 
elastic  parameters  for  these  units  have  been  adopted  from 
the  seismic  refraction  work  of  Shor  (1962)  and  Hales  and 
Asada  (1966)  and  the  microaftershock  array  work  of  Matumoto 
and  Page  (1969).  The  seismic  velocities  given  in  these 
studies  are  essentially  those  of  typical  crustal  and  upper 
mantle  material.  The  velocities  and  elastic  parameters  for 
these  units  are  listed  in  Table  3 •  1  •  Superposed  upon 
Figure  3-3  is  the  finite-element  grid  used  in  modeling  the 
fault  and  accompanying  dislocations.  The  grid  represents 
an  area  that  is  800  km  long  and  300  km  thick.  The  figure 
shows  the  Pacific  oceanic  plate  underthrusting  the  con¬ 
tinental  margin  beneath  the  eastern  Aleutian  arc.  The 
majority  of  the  material  modeled  in  this  finite  element 
grid  is  that  corresponding  to  the  oceanic  upper  mantle. 
Overlying  the  oceanic  upper  mantle  is  a  5  km  thick  zone  of 
oceanic  crust  which  also  underthrusts  the  continent  down 
to  a  depth  of  about  44  km.  Just  under  the  Alaska  trench 
we  have  inserted  a  thin  layer  of  typical  oceanic  sediments. 
The  fault  model  which  we  have  assumed  is  at  the  contact 
between  the  oceanic  crust  and  the  continental  crust.  The 
fault  starts  under  the  trench  with  a  dip  of  about  6°  and 
slowly  increases  its  dip  until  at  a  depth  of  28  km  the 
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fault  is  dipping  at  12°.  The  dip  continues  to  increase  so 
that  the  dip  is  15°  at  the  hypocenter  and  reaches  a  maximum 
of  20°  below  the  hypocenter. 

There  are  only  two  surface  faults  associated  with  this 
earthquake,  both  of  which  are  exposed  on  Montague  Island  — 
the  Patton  Bay  fault  and  the  Hanning  Bay  fault.  Geologic 
relations  (Plafker,  1967)  indicate  that  these  faults  are 
not  major  geologic  boundaries  but  rather  they  are  subsidiary 
to  the  zone  on  which  the  primary  faulting  motion  took  place. 
These  faults  can  be  considered  as  minor  imbrications  of  the 
megathrust.  Both  of  these  faults  have  been  mapped  to  strike 
approximately  parallel  to  the  continental  margin  and  the 
fault  motion  is  reverse  thrust  dipping  fairly  steeply  to 
the  northwest.  The  Patton  Bay  fault  has  a  large  component 
of  dip-slip  motion  associated  with  its  entire  length,  which 
extends  for  possibly  as  much  as  450  km  to  the  southwest 
(Plafker,  1972;  Malloy,  1964,  1965).  Reimnitz  (1966)  has 
inferred  that  this  fault  zone  extends  to  the  northeast  of 
Montague  Island  to  at  least  Hichinbrook  Island  some  50  km 
away.  The  strike-slip  component  is  measured  as  being  less 
than  one  meter  on  this  fault  so  that  the  motion  is  almost 
totally  dip-slip.  Von  Huene  et  al.  (1967)  carried  out 
seismic  and  echo  sounder  profiles  in  this  area  between 
Montague  Island  and  Kodiak  Island.  Their  results  indicate 
a  long  narrow  zone  of  faulting  with  the  vertical  attitude 
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of  the  fault  plane  estimated  to  be  60°.  By  the  observed 
deformation  of  the  sea  floor,  they  conclude  that  the  motion 
was  reverse  slip  along  this  steeply  dipping  plane  that  is 
inclined  landward.  This  fault  is  included  into  the  struc¬ 
tural  finite-element  model  as  a  reverse  fault  dipping  at 
58°  toward  the  continent.  This  fault  terminates  at  depth 
where  it  intersects  the  main  thrust  fault  at  a  depth  of 
about  25  km.  The  second  subsidiary  reverse  fault  observed 
on  Montague  Island,  the  Hanning  Bay  fault,  was  not  modeled 
in  this  study  because  of  the  short  length  (6  km)  of  the 
fault.  Another  high-angle  imbricate  reverse  fault  has  been 
proposed  to  break  the  surface  between  the  Patton  Bay  fault 
and  the  Aleutian  trench.  This  fault  has  been  inferred  to 
explain  the  large  vertical  displacements  on  Middleton 
Island.  However,  no  direct  physical  evidence  confirms  the 
existence  of  such  a  fault,  and  it  is  not  included  into  our 
model.  In  all,  a  total  of  26  nodes  in  the  finite-element 
grid  were  used  to  represent  the  megathrust  and  the  Patton 
Bay  fault,  21  nodal  elements  for  the  megathrust  and  5  nodal 
elements  for  the  subsidiary  fault. 

3 . 3  Static  Data. 

As  mentioned  in  the  introduction  to  this  chapter,  the 
crustal  deformation  accompanying  this  earthquake  was  very 
extensive.  Plafker  (1969)  has  described  in  detail  the 
regional  vertical  and  horizontal  displacements.  (See 


Plates  1  and  2  in  that  paper  for  detailed  contour  maps  of 
the  ground  deformation  and  the  location  of  the  observation 
sites.)  The  vertical  displacements  were  based  on  a  variety 
of  methods  of  measurements,  some  of  which  would  be  reliable 
only  if  the  net  vertical  deformation  was  large,  as  is  the 
case  for  this  event.  The  great  majority  of  the  measurements 
involved  measurements  of  the  movements  of  the  shoreline 
which  meanders  throughout  the  area  of  maximum  deformation. 
These  measurements  include  changes  in  tide  gauge  levels, 
measuring  the  change  in  the  upper  limit  of  barnacle  growth, 
direct  shoreline  changes,  etc..  Taken  individually,  these 
measurements  cannot  be  given  much  reliability,  however, 
when  the  entire  mass  of  these  observations  is  considered, 
including  correlation  between  geodetically  determined 
changes  in  bench  mark  levels,  the  data  become  quite  in¬ 
formative.  Plafker  (1969)  discusses  the  acquisition  of 
this  data  and  the  associated  estimate  of  the  errors  in¬ 
volved. 

Although  the  vertical  displacements  measured  after 
this  earthquake  were  large,  the  horizontal  displacements 
appear  to  be  even  larger  (Whitten,  196^4  ,  1965).  Unfor¬ 
tunately,  horizontal  displacements  do  not  lend  themselves 
to  the  ease  of  facility  of  measurement  as  do  the  vertical 
displacements  for  this  case.  Parkin  (1966)  has  described 
the  re tri angulation  network  that  was  occupied  after  the 
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earthquake.  Horizontal  displacements,  generally  in  the 
direction  of  the  seaward  motion  of  the  continent,  of  up  to 
20  m  were  observed.  Pope  (1972)  used  this  data  to  compute 
the  components  of  strain  on  the  surface.  The  surveys  to 
determine  horizontal  movements  are  too  poorly  controlled 
and  too  easily  subject  to  bias  to  enable  a  detailed  quali¬ 
tative  inversion  of  the  strains. 

In  this  chapter  we  will  limit  our  inversion  data  set 
to  vertical  displacements  only.  The  reason  for  this  is 
that  we  consider  the  vertical  displacement  data  to  be  much 
more  accurate  than  other  features  of  the  tectonic  deforma¬ 
tion  such  as  horizontal  shortening,  horizontal  displacements, 
and  changes  in  the  local  gravity  field.  The  vertical  dis¬ 
placement  data  are  taken  from  Plafker  (1965,  1969).  Since 
with  this  finite-element  method  we  are  limited  to  plane- 
sv.rain  problems  we  will  have  to  limit  our  data  set  to  points 
that  define  a  profile  perpendicular  to  the  strike  of  the 
megathrust.  We  chose  our  displacement  profile  to  coincide 
with  profile  BB'  in  Plafker's  papers  (1965,  1969,  1972). 

Only  one  major  surface  fault  Intersects  this  profile,  the 
Patton  Bay  fault  on  Montague  Island.  By  choosing  our  cross 
section  near  the  center  of  the  large  area  of  deformation, 
the  vertical  displacements  ar?  due  almost  totally  to  dip- 
slip  motion  on  the  fault ,  thus  contamination  of  the  data 
set  due  to  contributions  from  any  strike-slip  motion  is 
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minimized.  By  choosing  the  profile  in  this  position,  any 
effects  due  to  the  finite  length  of  the  fault  are  also 
minimized.  As  noted  above,  only  slight  amounts  of  strike- 
slip  motion  was  observed  along  this  cross  section  with  most 
of  it  being  on  the  subsidiary  reverse  faults  found  on 
Montague  Island.  This  absence  of  large  strike-slip  motion 
over  long  lengths  of  the  fault  allov?s  accurate  plane- 
strain  modeling  of  the  motions  involved.  We  also  restricted 
the  data  set  to  those  vertical  displacements  that  could  be 
confidently  projected  onto  this  profile.  Figure  3.4  shows 
this  cross  section  and  the  positions  of  the  data  available 
for  projection  onto  this  profile.  The  maximum  distance 
away  from  the  profile  of  a  data  point  was  about  75  km,  but 
about  9C%  of  the  available  data  points  were  within  40  km 
of  the  profile.  A  total  of  47  vertical  displacement  data 
points  were  chosen  along  the  profile  which  is  defined  for 
400  km  from  Middleton  Island  to  75  km  northwest  of  Cook 
Inlet.  Many  more  observations  were  available  within  the 
40  km  swath  on  either  side  of  the  profile,  however,  only 
those  points  that  were  not  near  a  curve  in  the  contours  or 
crossed  a  contour  were  acceptable  to  be  projected  onto  the 
profile.  The  projection  was  done  parallel  to  the  contours 
as  defined  by  Plafker  (1969).  This  projection  was  very 
close  (within  10°  in  most  instances)  to  a  perpendicular 
projection  onto  the  profile,  so  that  the  relative  location 


Figure  3.1|.  Location  of  the  vertical  displacement 
data  points  relative  to  the  profile. 


TABLE  3*2 


Distance 

Vertical 

Displacement 

(km) 

(m) 

101.0 

3-36 

102.0 

3.40 

182.0 

4.56 

184.0 

4.72 

184.5 

7.32 

185-0 

10.88 

190.0 

9.16 

191.0 

7-92 

192.0 

7-30 

196.0 

5-48 

207.0 

3.22 

211.0 

2.68 

212.0 

2.56 

214.0 

2.36 

219-0 

2.32 

226.0 

1.88 

229-0 

1.84 

231.0 

1.72 

232.0 

1.68 

233-0 

1.56 

235.0 

1.48 

245.0 

0.44 

247.0 

0.50 

258.0 

-0. 30 

Distance 

Ve rtical 
Displacement 

(km) 

(m) 

259-0 

-0.28 

285.0 

-1.20 

289.0 

-1.28 

291.0 

-1.40 

298.0 

-1.54 

300.0 

-1.62 

304.0 

-1.62 

335.0 

-1.72 

339.0 

-1.80 

342.0 

-1.75 

346.0 

-1.70 

350.0 

-1.64 

351.0 

-1.59 

355-0 

-1.52 

359-0 

-1 . 40 

362.0 

-0.92 

373.0 

-0.92 

401.0 

-0.30 

409.0 

-0.24 

413-0 

-0.24 

418.0 

-0.22 

423-0 

-0.24 

454.0 

0.00 

486.0 

0.44 

Table  3.2.  Observed  vertical  displacement  data  along 


profile  BB ’ . 
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of  the  data  points  on  this  profile  can  be  considered 
accurate  to  within  about  5  km  in  the  most  unfavorable 
cases.  The  corresponding  values  of  the  individual  pro¬ 
jected  data  points  on  the  profile  are  given  in  Table  3.2. 
These  values  and  their  respective  location  along  the  pro¬ 
file  will  appear  in  several  later  figures  in  this  chapter. 

The  origin  of  the  profile  is  some  100  km  southeast  of 
Middleton  Island.  For  reference,  the  most  southeasterly 
data  point  on  Middleton  Island  is  101.0  km  from  the  origin, 
and  the  profile  crosses  the  Patton  Bay  fault  at  a  distance 
of  185.0  km  from  the  origin  (B * ) .  The  sources  of  the  in¬ 
dividual  data  points  and  their  associated  errors  are  dis¬ 
cussed  elsewhere  (Plafker,  1969).  In  general,  the  data  are 
accurate  to  within  +0.3  m,  and  this  value  was  taken  in  the 
inversion  calculations. 

3.4  Calculated  Dislocation  Model. 

The  media  response  matrix.  A,  discussed  in  the  pre¬ 
vious  chapter  was  calculated  by  the  finite-element  tech¬ 
nique  for  the  structural  model  shown  in  Figure  3*3.  In 
this  technique,  the  static  displacement  on  the  nodal  seg¬ 
ments  at  the  free  surface  are  linearly  related  to  offsets 
imposed  on  the  designated  fault  nodes.  The  displacement  at 
every  one  of  the  nodal  segments  on  the  free  surface  due  to 
a  unit  offset  (1  m)  on  a  specified  fault  node  was  calculated. 
This  was  then  repeated  for  each  of  the  nodes  describing  the 
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fault.  system.  However,  since  the  observed  vertical  dis¬ 
placement  data  was  only  near  16  of  the  nodal  segments  on 
the  surface,  the  response  matrix  was  limited  to  those 
nodes.  Thus  we  have  defined  the  problem  of  estimating  the 
static  dislocation  on  26  fault  nodes  given  the  permanent 
static  offset  of  16  nodal  segments  located  on  the  free 
surface.  This  is  precisely  the  type  of  problem  that  was 
discussed  in  Chapter  2  for  which  we  formulated  the  sto¬ 
chastic  inversion  scneme  to  solve . 

In  this  problem,  the  operator  A  is  a  M  x  N  matrix, 
where  is  the  displacement,  calculated  at  the  point  on 
the  surface  where  the  ith  data  point  is  taken  due  to  a  unit 
dislocation  of  the  Jth  nodal  segment  of  the  fault.  Here 
M=?6  and  N=l6  Based  on  experience  in  calculating  best 
model  estimates  by  equation  (2.27),  it  was  found  that  much 
smoother,  hence  longer  wavelength,  solutions  were  calcu¬ 
lated  if  the  starting  model  was  some  "distance"  in  the 
model  space  away  from  the  null  model.  Therefore  we  chose 
to  use  Stauder  and  Bollinger's  (1966)  estimate  of  the  fault 
dislocation  as  the  starting  point  for  our  inversion.  This 
starting  model  turned  out  to  be  a  good  chcice  because  the 
inversion  scheme  smoothly  and  quickly  iterated  convergingly 
to  a  final  "best  fit"  model.  Just  to  make  sure  that  the 
final  model  that  we  obtained  was  not  wholly  dependent  on 
the  starting  model  that  we  chose,  we  then  repeated  the 
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inversion  using  Hastie  and  Savage's  (1970)  fault  disloca¬ 
tion  estimate  as  the  starting  model.  The  results  were  very 
similar  to  that  obtained  before.  We  therefore  feel  that 
this  final  model  is  not  very  dependent  on  the  starting 
model . 

The  upper  part  of  Figure  3-5  shows  this  vertical  sur¬ 
face  displacement  data  plotted  in  profile  and  the  calcu¬ 
lated  displacement  at  the  surface  nodes  of  the  finite 
element  grid.  The  fit  to  the  observed  data  is  extremely 
good  with  the  calculated  surface  displacement  field  fitting 
the  observed  data  used  in  the  inversion  to  within  a  RMS 
residual  of  about  3-1/2  cm,  and  the  fit  to  all  the  points 
in  the  data  set  is  not  far  from  this  value.  For  accuracy, 
only  those  data  points  which  were  very  near  a  surface  node 
in  the  finite  element  grid  were  used.  Thus,  out  of  the  set 
of  '*6  data  points  along  the  profile,  only  16  points  could 
be  actually  used  in  the  inversion.  An  increase  in  the 
number  of  surface  nodes  in  the  finite  element  grid  would 
probably  not  add  to  the  resolvability  o^  accuracy  of  tne 
slip  model,  since  the  limitations  In  these  quantities  were 
the  lack  of  soatial  coverage  of  the  data,  not  the  lack  of 
data  used.  The  slip  model  from  the  inversion  process  is 
shown  In  the  lower  half  of  Figure  3*5*  The  maximum  slip 
along  the  fault  is  33  m  at  a  point  below  Montague  Island. 

A  displacement  of  about  30  m  is  maintained  over  a  fault 
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width  of  about  60  km,  then  decreases  almost  linearly  at  a 
rate  of  0.3  m/km  over  the  next  100  km  as  the  fault  depth 
increases.  At  more  shallow  depths,  there  is  a  plateau  in 
slip  of  about  17  m,  which  would  correspond  to  the  fault 
surface  between  Middleton  Island  and  Montague  Island. 
However,  the  two  data  points  on  Middleton  Island  are  very 
important  in  this  model  in  that  their  values  almost  com¬ 
pletely  determine  the  amount  of  slip  along  the  top  150  km 
of  the  fault.  The  resolvability  of  this  plateau  will  be 
discussed  below.  The  slip  on  the  secondary  fault  is  not 
shown  in  this  figure,  but  it  averages  4  m  over  its  entire 
width  with  the  static  offset  on  the  node  at  the  surface 
constrained  to  be  equal  to  that  measured  for  the  scarp  on 
the  Patton  Bay  Fault  as  reported  by  Plafker  (1967).  The 
fault  offset  profile  on  the  main  fault  is  similar  in  shape 
to  that  proposed  by  Stauder  and  Bollinger  (1966)  who  used 
a  much  simpler  fault  model  and  ignored  the  effects  of 
geology . 

Integrating  the  area  under  the  slip  versus  fault 
width  curve,  we  find  that  we  have  an  average  slip  of  18.5 
m  over  a  260  km  fault  width.  This  slip  is  at  least  50? 
greater  than  that  predicted  by  Stauder  and  Bollinger  (1966), 
Savage  and  Hastie  (1966) ,  and  Hastie  and  Savage  (1970). 

One  check  to  see  if  the  average  dislocation  is  reasonable 
is  to  calculate  the  average  moment  and  compare  with  that 
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obtained  from  long  period  seismic  waves.  This  average 
moment  is  given  by 

M0  =  u  L  w  y 

where  u  is  the  average  fault  offset  (18.5  m) ,  L  the  lengtn 
(600  km),  W  the  width  (260  km),  and  y  is  the  average  rigid- 
ity  of  the  region  around  the  fault  (3*1  x  10  dyne/cm  ). 

By  using  the  rigidity  of  the  continental  crust,  the  mate¬ 
rial  in  which  most  of  the  deformation  takes  place,  we 

30 

obtain  an  average  moment  of  0.5  x  10  J  dyne-cm.  Kanamcri 
(1970)  arrives  at  a  moment  of  0.75  x  10^°  dyne-cm  on  the 
basis  of  long  period  (300  sec)  multiple  path  Love  and 
Rayleigh  waves.  At  these  long  periods,  the  surface  waves 
are  sampling  the  entire  fault  width  and  thus  should  give 
a  good  indication  of  the  average  moment.  These  two  values 
compare  very  favorably  indicating  that  indeed  there  were 
very  large  displacements  occurring  along  the  fault  sur¬ 
face.  B.  Minster  (personal  communication,  1973) »  on  the 
basis  of  a  systematic  inversion  of  world-wide  prate  motion 
data,  states  that  the  Pacific  plate  and  the  Alaskan  contin¬ 
ental  block  are  moving  relative  to  one  another  at  a  rate  of 
about  6  cm/year  at  the  location  of  our  profile. 

The  computed  average  slip  on  the  fault  leads  to  a  recur¬ 
rence  time  of  an  earthquake  of  this  magnitude  in  this  area 
of  once  every  300  years.  However,  if  the  central  portion 
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of  the  megathrust  with  its  average  30  m  slip  is  used  as 
representing  the  event ,  this  gives  a  recurrence  time  of 
500  years.  Pla^ker  and  Rubin  (1967)  obtain  a  repeat  time 
of  about  850  years  for  major  events  on  Middleton  Island 
based  on  the  radiometrically  determined  dates  of  a  set  of 
uplifted  marine  terraces  found  on  that  island.  However, 
Sykes  (1971)  has  expressed  great  uncertainty  about  estima¬ 
tions  of  recurrence  times  fox  major  events  in  this  region. 

Although  not  included  in  the  data  set  for  the  inver¬ 
sion,  the  measured  horizontal  displacement  field  was  ex¬ 
tensive.  Parkin  (1966)  gives  these  horizontal  movement 
vectors  which  are  made  with  a  free  adjustment  relative  to 
a  fixed  station  (Fishook  station)  located  about  14  km 
north  of  Palmer,  Alaska,  an  area  that  was  then  considered 
to  be  the  most  stable.  This  fixed  station  is  120  km  north¬ 
west  of  the  epicenter  of  the  main  shock.  As  in  the  case  of 
the  vertical  displacements,  only  those  horizontal  displace¬ 
ment  vectors  near  the  profile  line  were  chosen.  There  were 
23  of  these  vectors  in  the  vicinity  of  our  section.  These 
vectors  were  projected  onto  the  profile  and  their  component 
of  motion  in  the  direction  of  the  profile  taken.  The  re¬ 
sulting  displacements  are  shown  in  Figure  3*6.  The  da4- a 
points  nearest  the  fixed  station  are  the  most  accurate, 
being  first  order  surveys,  while  the  data  on  Montague  Island 
are  much  less  accurate,  being  based  on  third  order 


Figure  3.6.  Observed  and  calculated  relative  horizontal 
displacements  and  best  fit  dislocation  model.  A  positive 
horizontal  displacement  is  in  a  direction  E45°S. 
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observations  The  direction  of  motion  for  each  of  these 
points  shown  is  to  the  southeast.  ['he  horizontal  dis¬ 
placement  in  this  same  direction  is  calculated  from  the 
best  fit  slip  model  discussed  above  and  shown  in  the  figure. 
These  calculated  points  are  translated  relative  to  the  dis¬ 
placement  at  the  node  on  which  the  observed  apparent  zero 
isobase  is  projected.  The  resulting  displacements  form  a 
smooth  curve  except  for  the  irregularity  at  the  Patton  Bay 
Fault.  This  irregularity  is  not  resolvable  in  the  data 
shown  here.  Even  though  these  lateral  displacements  were 
not  used  in  the  inversion  scheme,  because  of  their  lack  of 
accuracy,  the  fit  is  surprisingly  good.  The  model  predicts 
a  movement  of  4  m  to  the  southeast  at  the  fixed  station.  A 
stable  area  for  displacement  reference  is  given  to  be  at 
least  120  km  farther  to  the  northwest  than  the  chosen  fixed 
station.  The  consistency  of  the  fit  to  both  the  horizontal 
and  vertical  displacement  data  seems  to  indicate  that  the 
model  geometry  that  was  initially  assumed  is  reasonably 
accurate . 

Figure  3.7  shows  a  contour  plot  of  the  calculated  dis¬ 
placement  field  in  two  dimensions  along  this  chosen  section. 
The  contour  values  are  indicated  on  the  figure  and  the  units 
are  in  km.  In  the  upper  half  of  this  figure  is  displayed 
the  calculated  two  dimensional  vertical  displacement  (Y 
direction  in  figure).  From  this  figure,  we  see  that  the 
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displacement  is  concentrated  under  Montague  Island  and 
above  the  fault  surface.  The  presence  of  the  Patton  Bay 
Fault  is  clearly  visible  on  this  plot.  This  partitioning 
of  the  displacement  field  is  due  to  the  effect  of  the 
nearby  free  surface.  In  the  lower  p^rt  of  this  figure,  we 
see  that  the  horizontal  displacement  (labeled  the  X  direc¬ 
tion  in  figure)  is  likewise  concentrated  immediately  above 
the  fault  surface. 

3.5  Resolvablllty  of  Features  in  the  Slip  Model. 

Since  the  data  used  in  the  inversion  are  not  perfectly 
accurate,  there  exist  model  perturbations  which  when  added 
to  our  best  fit  slip  model  would  still  fit  the  observed 
surface  displacement  data  to  some  chosen  degree  of  confi¬ 
dence.  If  we  can  estimate  the  errors  in  our  data,  then  we 
want  to  somehow  relate  these  errors  to  errors  in  cu1”  model. 
Such  a  relation  between  the  data  space  and  the  model  space 
exists  in  the  form  of  a  variance  operator  (equation  (2.33)) • 
This  operator  is  useful  in  this  application  in  the  follow¬ 
ing  manner.  If  we  take  some  perturbation,  6ra,  to  the 
calculated  slip,  then  this  perturbation  is  resolvable  by 
the  data  to  within  a  certain  confidence  interval  if  the 
following  inequality  holds, 

'Om*  V*  6m  >  k2(c)  . 
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V  Is  the  generalized  inverse  obtained  by  spectral  decom¬ 
position  of  the  variance  operator,  and  k(c)  is  the  coeffi¬ 
cient  associated  with  a  particular  confidence  interval. 

In  this  study  we  have  chosen  to  examine  model  perturbations 
at  the  95?  confidence  level,  so  that  in  this  case  the 
coefficient  associated  with  this  interval  is  1.96.  Using 
this  method,  we  can  test  chosen  perturbations  to  our  cal¬ 
culated  slip  model  and  compute  the  maximun  perturbation 
that  can  be  resolved  at  the  95?  confidence  level  by  the 
data.  We  note  that  these  tests  are  independent  of  the 
values  of  the  slip  model  itself,  and  only  pei-turbations  to 
this  model  can  be  checked  for  resolvability. 

The  variance  operator,  V,  for  this  case  is  a  26  x  26 
matrix.  The  generalized  inverse  of  this  matrix  is  found 
by  using  the  eigenvector  expansion  described  in  equation 
(2.32).  W  found  that  there  were  16  non-zero  eigenvalues 
associated  with  this  operator.  For  problems  where  the 
estimated  errors  are  very  small,  numerical  problems  may  be 
encountered  ir.  calculating  the  generalized  inverse  of  this 
operator.  These  numerical  problems  arise  from  the  fact 
that  round-off  errors  occur  in  the  computer  calculations 
of  the  eigenvalues.  For  small  eigenvalues,  the  problem  of 
distinguishing  non-zero  eigenvalues  from  the  zero  eigen¬ 
values  can  become  serious.  Fortunately,  this  is  not  the 
case  in  this  problem.  The  non-zero  eigenvalues  are  well 
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defined.  We  have  empirically  noted  that  the  number  of  non¬ 
zero  eigenvalues  of  the  variance  operator  V  is  equal  to 
the  number  of  independent  data  points  used  in  the  inver¬ 
sion.  The  inner  product  of  the  variance  operator  and  its 
generalized  inverse  form  a  projection  operator.  This  pro¬ 
jection  operator  is  then  checked  for  its  idempotent  proper¬ 
ties  to  make  sure  that  all  scaling  is  correct.  This  tost 
is  done  in  the  following  manner: 

(  VV+)(  VV+)  -  (  VV+)  <  E  (3.1) 

where  the  components  of  E,  ,  are  taken  to  be  some  small 
number  relative  to  the  :’ize  of  the  components  of  V. 

The  question  that  we  would  ultimately  like  to  answer 
with  a  study  of  this  type  is,  "What  is  the  maximum  pertu- 
bation  that  we  can  add  to  our  'best  fit*  slip  model  and 
still  satisfy  the  observed  data?"  Since  we  know  that  the 
size  of  the  maximum  perturbation  that  is  at  the  threshold 
of  detection  by  the  data  depends  on  the  distribution  of  the 
perturbation,  we  choose  three  perturbations  which  will 
elucidate  the  total  resolvability  of  our  slip  model.  We 
first  consider  how  much  of  a  slip  perturbation  we  can  add 
to  the  dislocations  in  the  hypocentral  region,  so  that  the 
rapid  fall-off  from  the  30  m  plateau  is  not  so  rapid. 

Figure  3-8a  shows  this  maximum  slip  perturbation.  The 


Figure  3.8.  Examples  of  the  resolvability  of  the  best  fit 
solution. 
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stippled  area  on  this  curve  is  the  maximum  slip  vhat  cculd 
be  added  in  this  region  and  still  be  undetected  by  the 
data.  It  is  seen  from  the  small  size  perturbation  in  this 
figure  that  the  data  control  very  closely  the  rate  of  fall- 
off  of  slip  in  thts  area  of  the  fault.  This  is  due  mainly 
to  the  fact  that  there  is  a  dense  network  of  data  points 
just  above  these  particular  nodes.  Next,  we  try  to  deter¬ 
mine  if  the  data  demand  the  existence  of  “he  17  m  plateau 
in  the  shallow  part  of  the  fault.  Figure  3 -8b  shows  the 
amount  of  slip  that  could  be  added  in  this  region.  We  see 
that  the  slip  gradient  in  this  region  could  be  smaller  than 
that  presented  in  our  best-fit  slip  model,  although  there 
still  appears  to  be  a  requirement  for  a  sharp  decay  in  slip 
up  the  fault  from  the  30  m  plateau.  The  slight  minimum  in 
slip  that  appears  in  this  region  o**  the  model  is  not  re¬ 
solvable  by  the  da^a.  In  Figure  3.8c  we  see  that  there  is 
almost  no  resolution  along  the  upper  part  of  the  fault. 

This  is  due  to  the  paucity  of  data  on  the  surface  above 
this  region.  In  order  to  explain  the  behavior  of  the  fault 
slip  in  this  region  we  have  to  appeal  to  arguments  based 
on  other  geophysical  data  than  the  statical  displacements. 
For  instance,  it  can  be  shown  that  large  fault  offsets  in 
the  area  of  the  trench  would  result  in  a  significant  amount 
of  strain  energy  stored  by  the  fault  in  that  region. 
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3 . 6  Averaging  Operators. 

In  Chapter  2  we  say  that  our  "best  fit"  estimate  of 
the  slip  model  is  in  reaUty  some  filtered  average  of  the 
true  slip  model.  TUs  filtering  operator  is  commonly  known 
as  the  averaging  operator.  Before  discussing  features  of 
our  final  model  it  is  to  our  advantage  to  know  the  extent 
of  the  avenging  that  is  taking  place  in  our  model.  The 
kernels  of  the  averaging  operator  are  taken  to  be  indivi¬ 
dual  rows  of  the  operator  matrix,  R,  as  defined  in  equation 
(2.29),  with  a  single  kernel  being  defined  for  each  fault 
element  comprising  the  total  fault  system.  We  note  here 
that  if  the  problem  is  linear,  as  it  is  in  this  case,  that 
these  kernels  do  not  depend  on'the  final  estimate  of  the 
"best  fit"  model. 

If  a  particular  slip  model  value  were  perfectly  well- 
known  by  the  inversion  then  the  averaging  component  cen¬ 
tered  on  that  fault  element  would  bQ  unity  and  all  the 
other  components  of  this  kernel  would  be  zero.  However,  in 
the  general  case  where  we  have  less  than  infinite  data  and 
the  data  that  we  do  have  are  somewhat  corrupted  by  noise, 
the  center  averaging  values  are  not  unity  and  the  other 
components  (off-diagonal  components  of  the  matiix  R)  are 
non-zero.  The  ability  to  resolve  the  details  of  the  actual 
dislocation  function  depends  on  two  features  of  this 
operator,  'tie  is  the  size  of  the  kernels.  This  depends 
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in  a.  general  way  on  the  availability  of  dcta  to  be  included 
into  the  inversion  that  are  sensitive  to  a  dislocation  ovei 
the  part  of  the  fault  model  that  we  are  testing.  As  the 
value  of  a  particular  diagonal  component  of  R  becomes  sub¬ 
stantially  less  than  unity,  our  ability  to  even  estimate 
the  slip  value  for  the  corresponding  model  component  de¬ 
creases.  The  other  facto’  is  the  averaging  width  of  the 
kernels.  This  averaging  width  is  expressed  by  tne  off- 
diagonal  elements  of  R.  If  these  off-diagonal  components, 
corresponding  to  the  fault  elemerts  "near"  the  particular 
11  element  we  are  examining  are  substantially  non-zero, 
ti.cn  the  estimated  "best  fit"  value  of  slip  that  we  obtain 
from  the  inversion  is  really  some  linearly  averaged  value 
of  the  actual  slip  values  in  the  vicinity  of  this  fault 
element.  These  ideas  are  probably  best  expressed  by 
examining  an  example  of  their  use. 

Figures  3.9  and  3.10  show  examples  of  the  averaging 
kernels  fov'  the  Alaska  earthquake  model.  The  coefficients 
of  the  rows  of  the  averaging  operator  are  shown  diagram- 
maticelly  at  the  position  of  the  respective  fault  node 
corresponding  to  the  components  of  this  row.  The  height 
of  the  bar  plotted  on  each  node  signifies  the  absolute 
value  of  the  averaging  coefficient  for  that  node.  For 
absolute  relerence,  in  Figure  3* 10c,  the  height  of  the 
outstanding  bar  is  0.997.  In  Figure  3-9,  we  have  plotted 


Figure  3.9.  Resolving  kernels  for  selected  nodal  segments 
along  the  megathrust.  View  is  a  perspective  of  the  mega- 
thrust  from  the  southeast  (left)  to  the  northwest  (light). 
Depths  and  profile  distances  are  in  km. 


Figure  3.10.  Resolving  kernels  for  selected  nodal 
segments  along  the  megathrust  (a»b^  and  subsidiary 
fault  (c,d).  View  is  a  perspective  of  the  megathrust 
from  the  southeast  (left)  to  the  ncrchwest  (right). 
Depths  and  profile  distances  are  in  km. 
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the  averaging  for  representative  nodes  along  the  megathrust 
presented  here  in  a  perspective  view.  The  arrow  in  the 
figure  indicates  on  which  node  the  averaging  is  centered. 

In  Figure  3.9a  we  see  that  for  the  upper  part  of  the  mega¬ 
thrust,  the  averaging  values  are  very  small,  again  showing 
our  lack  of  resolvability  in  this  area  of  the  fault.  In 
Figure  3 -9b  the  kernel  values  are  larger  in  amplitude  indi¬ 
cating  our  ability  io  estimate  the  slip  in  this  portion; 
however,  we  see  that  there  are  large  side  lobes.  The 
negative  averaging  coefficient,  indicated  by  the  bar  extend¬ 
ing  downward,  means  that  a  positive  dislocation  on  the 
centered  fault  node  could  be  traded  off  with  a  negative 
dislocation  on  this  node,  and  the  data  would  not  be  able  to 
tell  the  difference.  In  Figure  3 -9c  we  see  that  the  aver¬ 
aging  over  the  adjacent  nodes  to  either  side  of  the  cential 
node  is  fairly  severe,  and  that  there  is  a  slight  amount  of 
coupling  to  the  subsidiary  pault.  In  the  bottom  figure,  we 
see  that  the  amplitudes  start  to  become  more  peaked,  indi¬ 
cating  better  resolvability.  We  also  note  that  the  slip  In 
this  area  is  completely  uncoupled  from  the  clip  on  the  sub¬ 
sidiary  faulting.  This  shows  how  the  effect  of  the  subsid¬ 
iary  fault  is  very  localized  with  respect  to  the  megathrust. 

The  examples  are  continued  in  Figure  3.10.  Sections  a 
and  b  of  this  figure  continue  to  show  that  on  the  lower 
part  of  the  megathrust  we  are  ab±e  to  determine  fairly  well 
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the  "best  fit"  estimate  of  the  slip,  but  this  slip  is 
generally  averaged  over  the  one  or  two  adjacent  nodal  seg¬ 
ments.  Figure  3.10c  shows  the  averaging  kernel  for  the 
top  most  nodal  segment  of  the  representation  of  the  Patton 
Bay  fault  Here,  the  displacement  is  almost  exactly 
determined.  The  height  of  the  bar  is  almost  unity,  0.997, 
and  there  is  practically  no  spatial  averaging.  We  would 
expect  this  result,  considering  that  this  fault  segment 
breaks  the  surface  and  the  amount  of  dislocation  on  this 
nodal  segment  is  constrained  by  the  scarp  size  on  Montague 
Island.  Likewise  in  Figure  3-lOd,  the  dislocation  at  some 
depth  on  the  subsidiary  fault  ia  well  determined,  and  there 
is  practically  no  trade-off  in  dislocation  here  to  a  dis¬ 
location  on  the  megathrust. 

The  information  contained  in  the  averaging  operator 
can  be  summarized  by  defining  a  resolvabilxty  ratio  for 
each  kernel.  This  ratio  is  defined  as  the  ratio  of  the 
value  of  the  diagonal  coefficient  of  R  to  the  averaging 
half-width.  This  averaging  half-width,  though  somewhat 
ambiguous  in  some  Instances  because  or  asymmetries,  can 
usually  be  estimated, however .  The  averaging  half-width  is 
measured  from  the  central  fault  node  to  the  point  where  the 
averaging  first  crosses  zero.  This  ratio  is  convenient  and 
meaningful  in  the  sense  that  it  takes  into  account  both  the 
variables  involved  in  estimating  resolvability:  the  height 
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of  the  kernel  and  the  averaging  distance.  For  well  resolv¬ 
able  features  of  our  model  we  would  expect  large  ratios; 
for  less  well  resolvable  features,  smaller  ratios.  The 
results  for  this  particular  fault  model  are  shown  in  Figure 
3.11. 

In  this  figure,  we  see  that  for  the  upper  75  km  of  the 
megathrust,  there  is  a  total  lack  of  resolvability,  con¬ 
trolled  by  the  lack  of  data  which  are  sensitive  to  a  dis¬ 
location  in  that  area.  The  resolvability  is  slightly 
peaked  for  the  area  of  the  megathrust  immediately  under 
Middleton  Island,  but  again  there  is  no  resolvability  in 
the  area  between  the  islands.  The  resolvability  decreases 
rather  evenly  for  the  lower  end  of  the  megathrust.  This 
is  thought  to  be  due  to  the  fact  that  the  dislocations  are 
occurring  at  distances  farther  and  farther  away  from  the 
data,  thus  dislocation  averaging  starts  to  become  a 
problem  and  the  resolvability  is  reduced. 

3-7  Stress  and  Strain  Energy  Density  Change . 

In  terms  of  understanding  the  focal  processes  of  earth¬ 
quakes,  an  important  parameter  is  the  stress  drop.  In 
previous  studies  of  earthquake*,  the  stress  drops  were  ob¬ 
tained  through  empirical  formulas  or  exact  derivations  for 
special  purpose  geometry  of  the  crack  (for  example,  Starr, 
1928;  Knopoff,  3  958 i  Aki ,  1966).  The  stress  drop  over  some 
fault  dislocation  area  is  usually  given  by  the  following 
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type  relation, 


A  o 


nV 

W 


(3.2) 


where  y  is  the  rigidity,  Um  is  the  maximum  displacement,  W 
is  some  measure  of  the  size  of  the  fault,  and  n  is  some 
constant  dependent  upon  the  geometry  and  nature  of  the 
faulting.  Use  of  this  formula  results  in  stress  drop  values 
that  are  averaged  over  the  entire  fault  plane.  For  in¬ 
stance,  Brune  and  Allen  (1967)  estimate  the  stress  drop  from 
the  average  offset  given  by  Savage  and  Hastie's  (1966)  dis¬ 
location  of  the  1964  Alaska  earthquake  to  be  27  bars.  In 

this  calculation,  n  is  taken  to  be  1.33,  W=200  km,  Um= 

1 1  2 

13.3  m,  and  y=3.0  x  10"“  dyne/cm  .  If  we  were  to  use  their 
formulation  with  our  average  dislocation,  we  would  obtain 
a  stress  drop  of  30  bars.  Chinnery  (1969)  and  Sato  (1972) 
point  out  that  jn  order  to  evaluate  n  the  assumption  of  an 
infinite  length  fault  is  usually  made.  These  two  authors 
have  derived  the  expression  for  the  stress  drop  for  a 
finite  rectangular  fault,  and  they  show  that  the  stress 
drops  obtained  for  these  faults  are  smaller  than  what  one 
would  obtain  for  infinite  length  faults.  For  the  Alaska 
earthquake,  Sato  (1972)  estimates  n  to  be  0 . 9 7  when  a 
constant  displacement  over  the  fault  surface  in  an  ideal¬ 
ised  medium  is  considered.  His  resulting  estimate  of  the 


-95- 


stress  drop  using  this  parameter  and  our  average  disloca¬ 
tion  would  be  about  22  bars.  This  last  value  is  still  an 
estimate  of  the  stress  drop  averaged  over  the  entire  fault 
surface  .  However,  it  is  obvious  that  if  the  dislocation 
is  varying  over  the  fault  plane,  and  the  geometry  of  the 
pla^e  changes  with  distance,  the  stress  change  will  not  be 
a  constant  over  the  entire  fault  surface.  Jungel?  (1973) 
has  shown  that  for  several  earthquakes  the  stress  drop  can 
vary  along  the  fault  by  as  much  as  an  order  of  magnitude. 

To  estimate  the  stress  drop  along  the  width  of  the 
fault  we  can  apply  equation  (3*2)  to  each  of  the  nodal  seg 
merits  which  define  the  fault  plane,  but  there  is  always 
uncertainty  in  the  estimate  of  the  parameter  n.  Jungels 
(1973)  has  shown  a  more  direct  method  of  calculating  the 
fault  stress  drop  distribution  with  the  finite  element 

method.  This  is  accomplished  by  first  imposing  a  composit 

% 

pre3tress  field  on  the  structural  model.  This  is  done  by 
applying  a  dislocation  to  the  edges  of  the  structural 
model.  Prom  this  initial  state,  we  can  compute  the  equi¬ 
librium  final  state  that  would  be  caused  by  the  introduc¬ 
tion  of  our  best  fit  dislocation  mode] .  Then,  at  every 
point  of  the  structure,  the  difference  between  the  initial 
and  final  stress  fields  defines  the  stress  change.  If  A  7 
is  positive,  then  the  particular  model  caused  a  stress 
drop.  If  Ao  is  negative,  then  we  have  a  stress  increase. 
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it  is  clear  that  for  a  "dislocation"  model  the  magnitude  of 
the  stress  change  is  controlled  only  by  the  magnitude  of 
the  fault  offset  and  the  elastic  constants  of  the  struc¬ 
tural  model.  Thus  the  prestress  only  has  importance  in 
terms  of  the  strain  energy  change  where  the  sign  of  the 
stress  drop  matters.  The  shear  stress  change  approximately 
parallel  to  the  megathrust  is  calculated  by  this  method 
and  is  contoured  throughout  the  cross  section  studied. 
(Figure  3*12).  The  values  contoured  are  exact  away  from 
the  fault  plane,  but  for  the  nodes  defining  the  plane  it¬ 
self,  the  actual  stress  drop  is  approximately  twice  the 
value  shown.  This  error  arises  from  the  fact  that  a  linear 
behavior  of  displacement  is  assumed  inside  each  element  in 
the  finite-element  grid.  This  error  results  in  an  under¬ 
estimate  of  the  stress  change  on  the  fault  surface  ( Jungels , 
1973)-  Assuming  that  the  error  is  exactly  a  factor  of  2  in 
this  problem,  we  see  that  for  our  best  fit  model  the  stress 
change  along  the  fault  itself  varies  from  a  stress  increase 
of  86  bars  to  a  stress  drop  of  215  bars.  The  stress  change 
over  the  entire  width  of  the  fault  averages  a  stress  drop 
of  approximately  40  bars.  This  indicates  how  misleading  a 
value  of  the  average  stress  drop  could  be. 

The  details  of  the  stress  change  are  very  interesting. 
We  see  that  both  ends  of  the  fault  underwent  a  net  increase 
in  shear  stress.  For  the  shallow  portion  of  the  megathrust. 
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Figure  3.12.  Stress  change  and  strain  energy  released 
around  the  fault  surface  for  an  average  prestress 
level  equal  to  the  average  stress  drop.  Positive 
contour  values  of  stress  change  represent  stress 
drops.  Units  are  bars  and  10Z1  ergs  per  knw. 
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this  is  an  argument  in  favor  of  our  best  fit  solution  with 
its  small  offset  in  this  region.  It  follows  that  if  we 
increased  the  offset  on  the  shallow  end  of  the  fault,  even 
though  we  could  not  resolve  this  increase  with  the  surface 
data,  the  stress  change  would  increase  in  proportion  and 
this  in  turn  would  make  that  region  a  prime  candidate  for 
aftershock  activity.  The  fact  that  significant  aftershocks 
were  not  observed  here  (Algermissen  et  al . ,  1972)  argues 
for  our  best  fit  solution.  On  the  other  hand,  if  the  off¬ 
set  at  the  shallow  end  of  the  megathrust  were  the  maximum 
amount  indicated  on  Figure  3.8,  in  all  likelihood  the  dis¬ 
placements  would  rupture  the  free  surface.  If  this  was 
the  case,  the  stored  stress  would  be  relieved  and  thus 
there  would  he  very  little  or  no  aftershook  activity.  A 
search  of  the  literature  concerning  this  earthquake  re¬ 
vealed  that  there  seems  to  have  been  no  post-earthquaki 
reconnaissance  of  the  ocean  floor  in  the  vicinity  of 
Middleton  Island  and  further  toward  the  Aleutian  trench, 
so  that  the  possibility  of  this  occurring  cannot  be  ruled 
out.  A  hydrographic  and  ocean-bottom-scanning  sonar  survey 
of  the  area  to  the  southwest  of  Montague  Island  revealed 
fresh  scarps  on  older  en  echelon  faults  sub-parallel  to  the 
extension  of  the  Patton  Bay  fault  (Malloy  and  Merrill,  1969.) 
These  authors  attribute  these  scarps  to  the  Patton  Bay 
fault  system.  It  is  conceivable  that  much  of  the  strain 
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from  the  large  dislocations  on  the  very  shallow  end  of  the 
megathrust  could  be  relieved  in  this  fashion.  That  is, 
large  displacements  on  the  fault  surface  are  absorbed 
through  a  system  of  high  angle  bifurcations  of  the  main 
thrust  sheet.  We  will  see  in  the  next  chapter  that  this  is 
precisely  what  occurred  during  the  1971  San  Fernando, 
California  earthquake.  Unfortunately,  for  the  Alaska 
earthquake,  the  data  are  not  adequate  to  prove  or  disprove 
that  this  condition  existed,  and  further  speculation  along 
these  lines  seems  fruitless. 

Another  area  of  slight  stress  increase  on  this  figure 
is  found  in  that  region  where  the  offset  function  in  the 
best  fit  model  goes  through  a  local  minimum  between  the 
30  m  slip  plateau  and  the  17  m  slip  plateau.  We  have  seen 
from  the  above  discussion,  however,  that  this  minimum  is 
not  resolvable,  so  therefore  the  existence  of  the  stress 
Increase  in  this  region  is  not  resolvable  by  the  data. 

Most  of  the  stress  drop  along  the  fault  surface  occurs 
where  the  fault  dislocation  is  the  greatest.  The  maximum 
stress  drop,  215  bars,  is  found  along  the  megathrust  just 
below  the  intersection  of  the  Patton  Bay  Fault. 

A  plot  of  the  strain  energy  density  change  in  the 
media  as  seen  in  the  lower  half  of  Figure  3.12  illustrates 
that  most  of  the  energy  available  for  seismic  radiation 
would  come  from  the  central  area  of  the  megathrust  in  the 
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same  region  where  the  maximum  stress  drop  occurs.  The 

21  3 

contours  In  the  figure  are  In  units  of  10  ergs  Am  .  A 
direct  comparison  can  be  made  between  this  strain  energy 
density  plot  and  the  multiple  rupture  characteristics  that 
Wyss  and  Brune  ( 19 6 7 )  found  for  this  event.  These  authors 
interpret  the  P-wave  radiation  as  caused  by  a  multiple 
event  source  mechanism  whereby  the  rupture  initiating  at 
the  hypocenter  travels  up  the  fault  plane  triggering  dis¬ 
crete  seismic  events  larger  than  the  initial  event.  The 
largest  of  these  discrete  events  has  been  located  on  the 
megathrust  20  km  southeast  of  Montague  Island.  The  pulse 
from  this  region  was  delayed  from  the  initial  pulse  by  a 
time  corresponding  to  a  rupture  velocity  of  3.5  km/sec 
and  had  an  amplitude  significantly  larger  (up  to  30  times 
larger)  than  that  radiated  by  the  initial  snook.  This 
agrees  qualitatively  with  our  estimate  of  a  large  strain 
energy  density  change  of  up  to  0.28  x  10  ergs  Am  concen¬ 
trated  below  Montague  Island,  while  in  the  hypocentral 
region,  the  energy  density  change  is  computed  to  be  only 
0.02  x  10^°ergs/km^ . 

3.8  Accuracy  of  the  Plane-Strain  Approximation. 

We  would  like  to  somehow  approximate  the  errors  that 


occur  by  making  the  plane-strain  approximation  that  we  have 
taken  In  this  example.  One  way  of  getting  an  estimate  of 
this  error  is  to  approximate  the  fault  model  by  a  series  of 
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three- dimensional  Volterra  planes  and  to  use  the  static 
dislocation  theory  (Smylie  and  Mansinha,  1971;  or  equation 
(2.2))  for  a  three-dimensional  fault  in  a  homogenous  half¬ 
space.  Although  this  model  will  not  have  the  influences 
of  the  lateral  heterogenities  included,  it  will  serve  to 
estimate  how  good  or  bad  the  approximation  is  that  we  have 
made.  A  Volterra  approximation  to  the  finite-element 
structural  model  was  made.  This  model  consisted  of  22 
individual  fault  elements,  18  to  describe  the  megathrust 
end  4  to  describe  the  subsidiary  faulting.  The  Volterra 
fault  elements  are  planar  surfaces  centered  on  the  posi¬ 
tion  of  the  finite-element  fault  nodal  segments  and  extend¬ 
ing  halfway  to  the  adjacent  fault  nodal  segments.  Only 
those  fault  nodal  segments  were  modeled  on  which  there  was 
a  calculated  non-zero  displacement.  Several  of  the  nodal 
s&gments,  at  the  shallow  end  of  the  megathrust  and  at  the 
very  deep  end  of  the  megathrust,  had  "best  fit"  dislocation 
estimates  of  zero.  These  segments  were  not  modeled  with 
the  Volterra  approximations.  The  dislocation  which  is 
constant  over  the  planar  surfaces  was  taken  to  be  equal 
to  that  of  the  finite-element  fault  nodal  segment  at  the 
center.  The  parameters  for  this  model  approximation  are 
given  in  Table  3. 3* 

We  first  calculated  the  vertical  displacement  for  a 
profile  due  to  this  fault  model  with  the  length  of  each 


TABLE  3-3 


Fault  Segment  Dip  d  W  <Au> 


(deg) 

(km) 

(km) 

(m) 

1 

6.0 

8.5 

22.8 

0.99 

2 

6.0 

11.0 

21.4 

2.93 

3 

6.0 

13-0 

18.0 

10.17 

4 

6.0 

15 .0 

16.2 

16.72 

5 

6.0 

17.0 

18.2 

16.39 

6 

6.0 

19.0 

20.1 

14.66 

7 

6.0 

21.0 

16.0 

16.05 

8 

6.0 

22.5 

13.0 

24.91 

9 

6  .0 

24.0 

9-7 

29.98 

10 

6.0 

25.0 

7.1 

33-22 

11 

6.0 

26.0 

14.2 

29.42 

12 

6.0 

27-5 

20.0 

28.28 

13 

8.0 

29.5 

20.0 

29.99 

14 

11.5 

32.5 

20.4 

23.91 

15 

14.0 

36.5 

20.4 

15-35 

16 

15.0 

41.5 

20.8 

11.70 

17 

15.0 

47.0 

31.2 

2.17 

18 

15.0 

55.5 

41.6 

1.03 

19 

58.0 

0.0 

2.4 

4.05 

20 

58.0 

2.0 

4.0 

3.31 

21 

58.0 

5.5 

5-8 

7.18 

22 

58.0 

10.5 

8.0 

2.83 

Table  3.3.  Source  parameters  for  the  3-dimensional 
homogenous  approximation  to  the  finite-element  model 
of  the  Alaska  earthquake,  d  is  the  depth  to  the  top  of 
the  planar  fault  surface;  W  is  the  width  of  the  fault 
surface  measured  along  the  dip;  and  <Au>  is  the  fault 
dislocation. 
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planar  surface  taken  to  be  10,000  km,  or  effectively  in¬ 
finity,  to  thus  approximate  the  plane-strain  criterior.  The 
profile  was  taken  to  be  equi-distance  from  the  ends  of  the 
fault  elements  and  perp  .ndicular  to  the  strike  of  the 
system.  The  vertical  displacements  in  a  profile  were  then 
calculated  from  this  fault  system  but  now  the  lengths  of 
the  individual  fault  elements  were  set  to  600  km,  the 
approximace  lower  limit  for  the  fault  length  estimated  to 
be  appropriate  for  this  event.  The  profile  was  taken,  not 
across  the  center  of  the  fault  system,  but  at  a  position 
80  km  from  the  center  and  still  perpendicular  to  the  strike 
of  the  fault  system.  This  profile  is  220  km  from  one  end 
of  the  fault  and  380  km  from  the  other  end.  This  is  approx- 
mately  the  maximum  distance  profile  BB*  in  Figure  3-2  can 
be  considered  from  the  center  of  the  fault  system.  The 
estimated  errors  arising  from  the  plane-strain  approximation 
was  taken  to  be  the  difference  between  the  computed  dis¬ 
placements  for  these  two  profiles.  This  difference  is  a 
function  of  the  distance  away  from  the  origin  of  the  fault 
system.  The  origin  of  the  fault  system  is  taken  to  be  the 
point  at  which  the  shallow  end  of  the  megathrust  projects 
to  the  surface.  The  differences  are  presented  in  Figure 
3.13a.  It  is  seen  from  this  figure  that  the  maximum  dis¬ 
placement  error  expected  from  the  plane-strain  approxima¬ 
tion  would  be  about  0.35  m  for  this  particular  model.  The 


Figure  3*13-  Estimates  of  the  errors  to  the  cal¬ 
culated  vertical  displacements  (top)  and  horizon¬ 
tal  displacements  (bottom)  due  to  the  plane-strain 
assumption  as  a  function  of  distance  along  profile 
BB* . 
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differences  between  these  two  profiles  are  relatively  con¬ 
stant  at  about  0.25  m,  and  the  sign  of  the  error  is  such 
that  the  deformation  at  the  surface  is  being  underestimated. 
This  implies  that  the  free  surface  displacements  for  a 
fault  dislocation  model  with  a  finite  size  length  are 
slightly  larger  than  those  from  a  model  in  which  each  fault 
component  has  infinite  length.  Thus,  we  can  say  that  the 
displacements  calculated  for  the  finite-element  model  are 
an  upper  bound  to  that  necessary  to  fit  the  data.  Consid¬ 
ering  the  finiteness  of  the  length  of  the  actual  fault,  we 
would  need  only  slightly  less  displacement  on  the  mega¬ 
thrust. 

Now  the  horizontal  disp xacements  are  put  to  the  same 
test.  Horizontal  displacements  in  the  direction  perpendic¬ 
ular  to  the  strike  of  the  fault  system  from  the  Volterra 
dislocation  model  were  calculated  for  both  a  profile  due 
to  an  infinite  length  fault  and  for  a  profile  80  km  away 
from  the  center  of  a  fault  system  that  has  fault  element 
lengths  of  600  km.  Figure  3«13h  shows  the  differences 
between  the  former  and  the  latter  profiles.  It  is  seen 
here  that  as  the  profile  distance  becomes  greater  than 
half  the  fault  length,  the  errors  due  to  the  plane-strain 
approximation  start  to  become  more  significant.  It  is  seen 
here  that  for  the  farthest  distance  along  the  profile,  the 
expected  error  in  the  calculated  horizontal  displacements 
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ls  about  0.5  m.  The  sense  of  this  error  is  that  the  esti¬ 
mated  horizontal  displacements  made  under  the  plane-strain 
criterior  are  too  large.  Thus  in  Figure  3.6  where  we 
estimated  that  the  reference  station  for  the  measured  hori¬ 
zontal  displacements  (Fishook  station)  actually  moved  ^  m  to 
the  southeast,  we  have  to  revise  this  estimate  to  be  about 
3.5  m.  The  area  of  horizontal  stability;  that  is,  the  area 
where  no  horizontal  movement  was  expected,  is  still  some 
75-100  km  to  the  northwest  of  the  reference  station. 

We  will  now  briefly  examine  the  implications  of  the 
estimated  error  due  to  the  plane-strain  approximation. 

Since  only  the  vertical  displacements  were  actually  used 
in  the  inversion  procedure,  only  the  errors  associated  with 
these  measurements  will  affect  the  resolution  of  our  model. 
Since  the  estimated  errors  due  to  the  plane-strain  assump¬ 
tion  affecting  the  data  points  used  in  the  inversion  were 
about  equal  to  the  estimated  observational  error  of  the 
data  themselves,  we  can  estimate  that  at  most,  the  total 
variance  of  the  data  should  be  multiplied  by  a  factor  of  H. 
As  we  can  see  from  equation  (2.33),  if  we  want  to  recognize 
a  given  perturbation  to  our  model  at  the  same  confidence 
limit  as  before  (95?),  then  the  size  of  the  perturoation 
will  have  to  be  doubled.  This  means  that  in  Figure  3.8 
the  amplitude  of  the  stippled  area  will  be  doubled  if  we 
keep  the  shape  of  the  perturbation  as  before.  This  implies 
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that  the  steepness  of  the  rate  of  dislocation  fall-off  with 
distance  from  the  maximum  plateau  going  toward  the  hypo- 
center  is  not  quite  as  resolvable  as  before.  For  the  other 
two  perturbations  considered,  the  conclusions  arrived  at 
before  are  unchanged. 

3 . 9  Conclusions . 

A  dislocation  model  has  been  presented  for  the  1964 
Alaska  earthquake.  The  surface  displacements  from  this 
model  are  calculated  with  the  finite-element  numerical 
modeling  technique  in  which  the  effects  of  both  the  known 
geologic  heterogeneities  of  the  region  and  the  non-linearity 
of  the  assumed  fault  plane  are  taken  into  account.  The 
dislocation  model,  which  was  obtained  using  a  stochastic 
inversion  scheme,  fits  with  high  precision  both  the  ob¬ 
served  vertical  and  horizontal  displacements.  The  calcu¬ 
lated  static  offset  along  the  fault  plane  was  found  to  be 
variable  and  to  have  a  maximum  amplitude  much  greater  than 
previously  imagined,  although  the  average  moment  agrees 
with  that  observed  from  long  period  surface  waves.  The 
two-dimensional  displacement  field  was  found  to  be  strongly 
partitioned  above  and  below  the  fault  surface,  with  most  of 
the  displacement  occurring  above  the  fault.  The  calculated 
displacement  at  the  shallow  end  of  the  fault  model  was 
found  to  be  almost  non-resol vable  due  to  the  lack  of  sur¬ 
face  displacement  data,  while  the  displacement  near  the 


hypocenter  was  well  constrained  by  the  data.  Along  with 
the  displacement  calculated  along  the  fault  surface,  both 
the  stress  drop  and  the  strain  energy  density  varied  widely 
The  maximum  stress  drop  found  was  218  bars ,  while  at  both 
ends  of  the  fault  the  stress  field  increased  as  a  result  of 
the  static  dislocations.  The  region  of  maximum  stress  drop 
and  maximum  strain  energy  density  change  calculated  from 
this  static  study  was  found  tc  correspond  to  the  region  of 
maximum  compressional  wave  radiation.  The  errors  caused 
by  the  plane  strain  approximation  for  this  event  were 
analyzed  and  found  not  to  affect  any  of  the  above  conclu- 


